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Preface to Matrix Analysis 
Matrix Analysis 


Bellman has called matrix theory 'the arithmetic of higher mathematics.’ 
Under the influence of Bellman and Kalman, engineers and scientists have 
found in matrix theory a language for representing and analyzing 
multivariable systems. Our goal in these notes is to demonstrate the role of 
matrices in the modeling of physical systems and the power of matrix 
theory in the analysis and synthesis of such systems. 


Beginning with modeling of structures in static equilibrium we focus on the 
linear nature of the relationship between relevant state variables and express 
these relationships as simple matrix-vector products. For example, the 
voltage drops across the resistors in a network are linear combinations of 
the potentials at each end of each resistor. Similarly, the current through 
each resistor is assumed to be a linear function of the voltage drop across it. 


And, finally, at equilibrium, a linear combination (in minus out) of the 
currents must vanish at every node in the network. In short, the vector of 
currents is a linear transformation of the vector of voltage drops which is 
itself a linear transformation of the vector of potentials. A linear 
transformation of n numbers into m numbers is accomplished by 
multiplying the vector of n numbers by an m-by- n matrix. Once we have 
learned to spot the ubiquitous matrix-vector product we move on to the 
analysis of the resulting linear systems of equations. We accomplish this by 
stretching your knowledge of three-dimensional space. That is, we ask what 
does it mean that the m-by- n matrix X transforms {%” (real n-dimensional 
space) into SX? We shall visualize this transformation by splitting both i” 
and SR” each into two smaller spaces between which the given X behaves 
in very manageable ways. An understanding of this splitting of the ambient 
spaces into the so called four fundamental subspaces of X permits one to 
answer virtually every question that may arise in the study of structures in 
static equilibrium. 


In the second half of the notes we argue that matrix methods are equally 
effective in the modeling and analysis of dynamical systems. Although our 
modeling methodology adapts easily to dynamical problems we shall see, 
with respect to analysis, that rather than splitting the ambient spaces we 
shall be better served by splitting X itself. The process is analogous to 
decomposing a complicated signal into a sum of simple harmonics 
oscillating at the natural frequencies of the structure under investigation. 
For we shall see that (most) matrices may be written as weighted sums of 
matrices of very special type. The weights are eigenvalues, or natural 
frequencies, of the matrix while the component matrices are projections 
composed from simple products of eigenvectors. Our approach to the 
eigendecomposition of matrices requires a brief exposure to the beautiful 
field of Complex Variables. This foray has the added benefit of permitting 
us a more careful study of the Laplace Transform, another fundamental tool 
in the study of dynamical systems. 


--Steve Cox 


Matrix Methods for Electrical Systems 


Nerve Fibers and the Strang Quartet 


We wish to confirm, by example, the prefatory claim that matrix algebra is 
a useful means of organizing (stating and solving) multivariable problems. 
In our first such example we investigate the response of a nerve fiber to a 
constant current stimulus. Ideally, a nerve fiber is simply a cylinder of 
radius a and length / that conducts electricity both along its length and 
across its lateral membrane. Though we shall, in subsequent chapters, delve 
more deeply into the biophysics, here, in our first outing, we shall stick to 
its purely resistive properties. The latter are expressed via two quantities: 


1. p;, the resistivity in Qcm of the cytoplasm that fills the cell, and 
2. Pm, the resistivity in Qcm? of the cell's lateral membrane. 


A 3 compartment model of a nerve cell 


Although current surely varies from point to point along the fiber it is hoped 
that these variations are regular enough to be captured by a 
multicompartment model. By that we mean that we choose a number NV and 
divide the fiber into N segments each of length 4: Denoting a segment's 


axial resistance 


= Pin 
Tra 
and 
membrane resistance 
Pm 
y= a a 
27 an 


we alrive at the lumped circuit model of [link]. For a fiber in culture we 
may assume a constant extracellular potential, e.g., zero. We accomplish 
this by connecting and grounding the extracellular nodes, see [link]. 

A rudimentary circuit model 


i i RK, 
Dee ee 3 ' 


[link] also incorporates the exogenous disturbance, a current stimulus 
between ground and the left end of the fiber. Our immediate goal is to 
compute the resulting currents through each resistor and the potential at 
each of the nodes. Our long--range goal is to provide a modeling 
methodology that can be used across the engineering and science 
disciplines. As an aid to computing the desired quantities we give them 
names. With respect to [link], we label the vector of potentials 


PSH we we. Gy) 


and the vector of currents 


y=(y1 Yo Ys Ys Ys Yo). 


We have also (arbitrarily) assigned directions to the currents as a graphical 
aid in the consistent application of the basic circuit laws. 
The fully dressed circuit model 


We incorporate the circuit laws in a modeling methodology that takes the 
form of a Strang Quartet: 


e (S1) Express the voltage drops viae = — (Aza). 

e (S2) Express Ohm's Law via y = Ge. 

e (S3) Express Kirchhoff's Current Law via A’ y = —f. 
e (S4) Combine the above into A7GAz = f. 


The A in (S1) is the node-edge adjacency matrix -- it encodes the 
network's connectivity. The G in (S2) is the diagonal matrix of edge 
conductances -- it encodes the physics of the network. The f in (S3) is the 
vector of current sources -- it encodes the network's stimuli. The 
culminating A?GA in (S4) is the symmetric matrix whose inverse, when 
applied to f, reveals the vector of potentials, az. In order to make these 
ideas our own we must work many, many examples. 


Example 


Strang Quartet, Step 1 


With respect to the circuit of [link], in accordance with step (S1), we 
express the six potential differences (always tail minus head) 


€1 — %1— 
€2 — 12 
€3 —= %2— #3 
€4 = 13 
€5 = @3— X4 
€6 — L4 


Such long, tedious lists cry out for matrix representation, to wit 
e = — (Az) where 


=—l1 1 0: 
U.. = 0 
Re 0 -1 1 0 
0 0 =-1 OQ 
0 O -1 1 
CY @ =! 


Strang Quartet, Step 2 


Step (S2), Ohm's Law, states: 
Ohm's Law 


The current along an edge is equal to the potential drop across the edge 
divided by the resistance of the edge. 
In our case, 


ej, ej, 
Up Tyas? and ae Ae i 
or, in matrix notation, y = Ge where 
1 
BR 0 0 0 0 0 
1 
0 Re 0 0 0 0 
1 
An 0 0 RB 0 0 0 
— 1 
0 0 0 Re 0 0 
0 0 0 0 & 0 
0 0 0 0 0 3 


Strang Quartet, Step 3 


Step (S3), Kirchhoff's Current Law, states: 


Kirchhoff's Current Law 


The sum of the currents into each node must be zero. 
In our case 


m—-y =O 
Yi — y2— 3 = 0 
¥3 — Ya — Ys = O 

Ys — ¥6 = 0 


or, in matrix terms 


where 


—1 0 0 0 0 90 10 


— 
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Strang Quartet, Step 4 


Looking back at A: 
—-1 1 0 0 
0 -1 0 
0 -1 1 O 
A= 
0 O -1 O 
0 O -1 1 
0 0 O -1 


we recognize in B the transpose of A. Calling it such, we recall our main 
steps 


- (Sle = — (Az), 
e (S2) y = Ge, and 
¢ (S3) Aly = —f. 


On substitution of the first two into the third we arrive, in accordance with 
(S4), at 
Equation: 


ATG Az = f. 


This is a system of four equations for the 4 unknown potentials, x; through 
x4. As you know, the system [link] may have either 1, 0, or infinitely many 
solutions, depending on f and A?G.A. We shall devote (FLX ME CNXN 
TO CHAPTER 3 AND 4) to an unraveling of the previous sentence. For 


now, we cross our fingers and ‘solve’ by invoking the Matlab program, 
fibl.m . 


Results of a 64 compartment simulation 
Potential along 64 compartment fiber 
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Results of a 64 compartment simulation 
Axial current along 64 compartment fiber 


Membrane current along 64 compartment fiber 
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This program is a bit more ambitious than the above in that it allows us to 
specify the number of compartments and that rather than just spewing the x 
and y values it plots them as a function of distance along the fiber. We note 
that, as expected, everything tapers off with distance from the source and 
that the axial current is significantly greater than the membrane, or leakage, 
current. 


Example 


We have seen in the previous example how a current source may produce a 
potential difference across a cell's membrane. We note that, even in the 
absence of electrical stimuli, there is always a difference in potential 
between the inside and outside of a living cell. In fact, this difference is the 
biologist's definition of ‘living.’ Life is maintained by the fact that the cell's 
interior is rich in potassium ions, K*, and poor in sodium ions, Na*, while 
in the exterior medium it is just the opposite. These concentration 
differences beget potential differences under the guise of the Nernst 
potentials: 


Nernst potentials 


RT [Nal], RT [K], 
— log and Kx = —— log 
F [Na], F [K], 


u 4 


where £ is the gas constant, T’ is temperature, and F’ is the Faraday 
constant. 


Associated with these potentials are membrane resistances 
Pm,Na and Pm,K 


that together produce the p,, above via 


1 1 | 
ss ee be 
Pm Pm,Na Pm,K 


and produce the aforementioned rest potential 


Ena Zz Ex ) 
Pm,Na Pm,Na 


Bm = Pm ( 


With respect to our old circuit model, each compartment now sports a 
battery in series with its membrane resistance, as shown in [link]. 
Circuit model with resting potentials 


Revisiting steps (S1-4) we note that in (S1) the even numbered voltage 
drops are now 

€9 = 209 — Eq 

€e4=—23— Ey 

€g = 414 — En 
We accommodate such things by generalizing (S1) to: 


e (S1') Express the voltage drops as e = b — Az where 6 is the vector 
of batteries. 


No changes are necessary for (S2) and (S3). The final step now reads, 


e (S4') Combine 


Returning to [link], we note that 


Em 


A'Gb= — 
and Gb R,. 


ee Oo - CO KF © 
Se eS © 


This requires only minor changes to our old code. The new program is 
called fib2.m and results of its use are indicated in the next two figures. 


Results of a 64 compartment simulation with batteries 
Potential along 64 compartment fiber 
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Results of a 64 compartment simulation with batteries 
Axial current along 64 compartment fiber Membrane current along 64 compartment fiber 
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Exercises for Matrix Methods for Electrical Systems 
Exercise: 


Question 1 


In order to refresh your matrix-vector multiply skills please calculate, by 
hand, the product A?G‘A in the 3 compartment case and write out the 4 
equations in the vector equation we arrived at in step (S4): ATGAa = f. 


Feedback 
The second equation should read 
Equation: 
—2%, + 229 — 23 r2 =) 
R; Re. 
Exercise: 
Question 2 


We began our discussion with the ‘hope’ that a multicompartment model 
could indeed adequately capture the fiber's true potential and current 
profiles. In order to check this one should run fib1.m with increasing values 
of N until one can no longer detect changes in the computed potentials. 


e (a) Please run fibi.m with N = 8, 16, 32, and 64. Plot all of the 
potentials on the same (use hold) graph, using different line types for 
each. (You may wish to alter f1b1.m so that it accepts N as an 
argument). 


Let us now interpret this convergence. The main observation is that the 
difference equation, [link], approaches a differential equation. We can see 
this by noting that 


u(z) = ~ 


acts as a spatial 'step' size and that xx, the potential at (k — 1)d/(z), is 
approximately the value of the true potential at (k — 1)d/(z). Ina slight 
abuse of notation, we denote the latter 


x((k — 1)dl(z)) 


Applying these conventions to [link] and recalling the definitions of R; and 
R,, we see [link] become 


2 
a” —2x(0) + 2x(d(z)) — x(2d(z 2radl(z 
a? ~2(0) + 20(al(2)) ~ 2(2d(2)) | 2radl2) gy.) _ 4 

Pi al(z) Pm 
or, after multiplying through by qed) , 

apm —x(0) + 2x(A(z)) — x(2d(z)) 
Pi al(z*) 

. We note that a similar equation holds at each node (save the ends) and that 
as NV — oo and therefore a/(z) — 0 we arrive at 


+ 2x(di(z)) =0 


Equation: 
d? 20; 
——2(z) — th t(z)=0 
d 22 aPm 
e (b) With wp = ma ~ show that 
Equation: 


(2) = asinh(/2y2) + B cosh ( 22) 


satisfies [link] regardless of a and £3. 


We shall determine a and ( by paying attention to the ends of the fiber. At 
the near end we find 


which, as Q(z) — 0 becomes 
Equation: 


_ ito 


an Ta? 


a 
dz 


At the far end, we interpret the condition that no axial current may leave the 
last node to mean 
Equation: 


qet) —0 


e (c) Substitute [link] into [link] and [link] and solve for a and @ and 
write out the final x(z). 

e (d) Substitute into x the /, a, p;, and p,, values used in fib1.m, plot the 
resulting function (using, e.g., eZpLOt) and compare this to the plot 
achieved in part (a), 


Matrix Methods for Mechanical Systems: A Uniaxial Truss 


Introduction 


We now investigate the mechanical prospection of tissue, an application 
extending techniques developed in the electrical analysis of a nerve cell. In 
this application, one applies traction to the edges of a square sample of 
planar tissue and seeks to identify, from measurement of the resulting 
deformation, regions of increased hard ¥ ss/ or stiffness.' For a sketch of 
the associated apparatus, visit the Biaxial Test site . 


A Uniaxial Truss 


A uniaxial truss 


As a precursor to the biaxial problem let us first consider the uniaxial case. 
We connect 3 masses with four springs between two immobile walls, apply 
forces at the masses, and measure the associated displacement. More 
precisely, we suppose that a horizontal force, f;, is applied to each m,, and 
produces a displacement x;, with the sign convention that rightward means 
positive. The bars at the ends of the figure indicate rigid supports incapable 
of movement. The k; denote the respective spring stiffnesses. The analog of 
potential difference (see the electrical model) is here elongation. If e; 
denotes the elongation of the 7th spring then naturally, 


Cp) — £] 
€2 — %.—- TZ] 


€3 — %3—- 2 


€4 = —23 


or, in matrix terms, e — Aa, where 


1 O 0 

—— £ 
A= 

y =) 2 

0 O -1 


We note that e; is positive when the spring is stretched and negative when 
compressed. This observation, Hooke's Law, is the analog of Ohm's Law in 
the electrical model. 


Hooke's Law 
The restoring force in a spring is proportional to its elongation. We call 
the constant of proportionality the stiffness, k;, of the spring, and 
denote the restoring force by yj. 
The mathematical expression of this statement is: y; = kje,, or, 
in matrix terms: y = Ke where 


k 0 0 0 
xra|0 m0 0 
0 0 ks 0 
0 0 0 kg 


The analog of Kirchhoff's Current Law is here typically called “force 
balance.’ 


force balance 
Equilibrium is synonymous with the fact that the net force acting on 
each mass must vanish. 
In symbols, 


Yi-yo—fi =0 
y2 — y3 — fe =0 


y3 — ys — fg = 0 


or, in matrix terms, By = f where 


fi 1-1 0 0 
f= |felandB=|0 1 -1 O 
fs 0 0 1 -!1 


As in the electrical example we recognize in B the transpose of A. 
Gathering our three important steps: 


Equation: 

e= Ax 
Equation: 

y= Ke 
Equation: 

A'y=f 


we alrive, via direct substitution, at an equation for x. Namely, 
(A*y = f) > (A’ Ke = f) = (A’ KAz = f) 


Assembling A? K Az we arrive at the final system: 
Equation: 


kitk. -—ke 0 Ly fi 
—kp kot+k 3 —kg ro) = | fe 
0 —k3 kg+ka] \x3 fr 


Gaussian Elimination and the Uniaxial Truss 


Although Matlab solves systems like the one above with ease our aim here 
is to develop a deeper understanding of Gaussian Elimination and so we 
proceed by hand. This aim is motivated by a number of important 
considerations. First, not all linear systems have solutions and even those 
that do do not necessarily possess unique solutions. A careful look at 
Gaussian Elimination will provide the general framework for not only 
classifying those systems that possess unique solutions but also for 
providing detailed diagnoses of those defective systems that lack solutions 
or possess too many. 


In Gaussian Elimination one first uses linear combinations of preceding 
rows to eliminate nonzeros below the main diagonal and then solves the 
resulting triangular system via back-substitution. To firm up our 
understanding let us take up the case where each k; = 1 and so [Link] takes 
the form 


Equation: 
2 —]l 0 “1 fi 
—l 2 —] “v2 = fo 
0 —1 2 L3 fs 


We eliminate the (2,1) (row 2, column 1) element by implementing 


1 
new row 2 = old row 2 + 73 Tow 1 


bringing 
2-1 0\ /2, fi 
e -l||e2|=l|A+ fh 
0 -1 2 £3 fs 


We eliminate the current (3,2) element by implementing 


2 
new row 3 = old row 3 + re 2 


bringing the upper-triangular system 


2-1 0\ fz, fi 
4 2 
0 0 4) \zs fg+72 +4 


One now simply reads off 


fi + 2fo+ 3fs 


L3 = A 


This in turn permits the solution of the second equation 


2(w+ht+F)  f+oh+ fs 


VES 3 9 


and, in turn, 


_#ath _ sitthth 


2 4 


One must say that Gaussian Elimination has succeeded here. For, regardless 
of the actual elements of f, we have produced an 2 for which 
A’ K Az = f. 


Alternate Paths to a Solution 


Although Gaussian Elimination remains the most efficient means for 
solving systems of the form Sa = f it pays, at times, to consider alternate 
means. At the algebraic level, suppose that there exists a matrix that 
‘undoes' multiplication by S in the sense that multiplication by 2~' undoes 
multiplication by 2. The matrix analog of 2-12 = 1 is 


S'ts=lI 


where J denotes the identity matrix (all zeros except the ones on the 
diagonal). We refer to S~? as: 


Inverse of S 
Also dubbed "S inverse" for short, the value of this matrix stems from 
watching what happens when it is applied to each side of Sa = f. 
Namely, 


Se=f)=(S Se={S"f)=S Uz2=S "f) > (e=S fF) 


Hence, to solve Sx = f for z it suffices to multiply f by the inverse of S. 


Gauss-Jordan Method: Computing the Inverse of a Matrix 


Let us now consider how one goes about computing S'—!. In general this 
takes a little more than twice the work of Gaussian Elimination, for we 
interpret 


Sts=I 


as n (the size of S) applications of Gaussian elimination, with f running 
through n columns of the identity matrix. The bundling of these n 
applications into one is known as the Gauss-Jordan method. Let us 
demonstrate it on the S appearing in [link]. We first augment S with J. 


2-1 0/100 
-1 2 -1 | 010 
001 


We then eliminate down, being careful to address each of the three f 
vectors. This produces 


2-10 | 1 00 
0 2 -1 | 410 
00 3 lai 


Now, rather than simple back--substitution we instead eliminate up. 
Eliminating first the (2,3) element we find 


221. | L @ 0 
Org OL g a a 
oo 7/331 
Now, eliminating the (1,2) element we achieve 
200/21 4 
Og Be fig ae 
003/351 


In the final step we scale each row in order that the matrix on the left takes 
on the form of the identity. This requires that we multiply row 1 by 4, TOW 


2 by 3, and row 3 by 3, with the result 


Co CO 
oOo -_- © 
Ke CO © 
AlR welR Aloo 
NO ) 
Blo we Ale 


Now in this transformation of S into J we have, ipso facto, transformed I to 
S—1; i.e., the matrix that appears on the right after applying the method of 
Gauss-Jordan is the inverse of the matrix that began on the left. In this case, 


ae 


Ale wle Ble 
mle RR wl 
mle mle Ale 


One should check that S~! f indeed coincides with the 2 computed above. 


Invertibility 


Not all matrices possess inverses: 


singular matrix 
A matrix that does not have an inverse. 


Example: 
A simple example is: 


oy 


Alternately, there are 


Invertible, or Nonsingular Matrices 
Matrices that do have an inverse. 


Example: 
The matrix S that we just studied is invertible. Another simple example is 


i 4) 


Matrix Methods for Mechanical Systems: A Small Planar Truss 


We return once again to the biaxial testing problem, introduced in the 
uniaxial truss module. It turns out that singular matrices are typical in the 
biaxial testing problem. As our initial step into the world of such planar 
structures let us consider the simple truss in the figure of a simple swing. 
A simple swing 


We denote by 2, and £2 the respective horizontal and vertical displacements 
of m, (positive is right and down). Similarly, f; and f2 will denote the 
associated components of force. The corresponding displacements and 
forces at mg will be denoted by x3, x4 and fs, f4. In computing the 
elongations of the three springs we shall make reference to their unstretched 
lengths, Ly, Le, and Lz. 


Now, if spring 1 connects {0, —L;} to {0, 1} when at rest and {0, —L,} to 


{x1, £2} when stretched then its elongation is simply 
Equation: 


ej; = [2012 (xo + In)’ = Ly 


The price one pays for moving to higher dimensions is that lengths are now 
expressed in terms of square roots. The upshot is that the elongations are 


not linear combinations of the end displacements as they were in the 
uniaxial case. If we presume, however, that the loads and stiffnesses are 
matched in the sense that the displacements are small compared with the 
original lengths, then we may effectively ignore the nonlinear contribution 
in [link]. In order to make this precise we need only recall the 

Taylor development of the square root of (1 + t) 


The Taylor development of 21 + t¢ about t = 0 is 
—— t 
v2i+t=1+ 5 +0(#) 


where the latter term signifies the remainder. 
With regard to e; this allows 


Equation: 
e, = V 2a? +a? + 2aL,+L2-L, 
= [y4/21+ ae + oe — In 
Equation: 


2 
e, = E,+ ons ++ Bo ++ 1,0( (3 + 221 | ) — fy 


xy+297 xy +297 2x2 A 
= 2557 = Bia + L,O “Te al T. 


If we now assume that 
Equation: 


2 
xy + £2 


is small compared to x 
Ly p 2 


then, as the O term is even smaller, we may neglect all but the first terms in 
the above and so arrive at 


€1 = r2 


To take a concrete example, if £; is one meter and x, and 22 are each one 
2 2 
Ly +22 


centimeter, then x2 is one hundred times —5 Li 


With regard to the second spring, arguing as above, its elongation is 
(approximately) its stretch along its initial direction. As its initial direction 
is horizontal, its elongation is just the difference of the respective horizontal 
end displacements, namely, 


€2 — %3 —- @ 


Finally, the elongation of the third spring is (approximately) the difference 
of its respective vertical end displacements, i.e., 


€3 = @4 
We encode these three elongations in 


0 1 
e=— Ax where A=]-1 0 
0 O 


oF © 


0 
0 
1 


Hooke's law is an elemental piece of physics and is not perturbed by our 
leap from uniaxial to biaxial structures. The upshot is that the restoring 
force in each spring is still proportional to its elongation, i.e., y; = k,e; 
where k, is the stiffness of the jth spring. In matrix terms, 


ky O O 
y= Ke where K=]|0 k O 
0 O kg 


Balancing horizontal and vertical forces at m , brings 
(—ye) — fi =0 


and 


w— f2=0 


while balancing horizontal and vertical forces at mz brings 


y2— fe =0 
and 
y3 — fa =0 
We assemble these into 
G. =1. 0 
1 0 O 
By=f where B= 0 1 ol’ 
0 0 1 


and recognize, as expected, that B is nothing more than Af. Putting the 
pieces together, we find that a must satisfy Sw = f where 


ko 0 —ko O 

0 k O O 
S-A'KA= 

—ko 0 ky O 

0 0 O ks 


Applying one step of Gaussian Elimination brings 


ko O —ke O Ly fi 

Ok 0 O|}fm) | fh 
0 0 O}fa3} | fit fs 
0 0 kg “4 fa 


and back substitution delivers 


0=fit fs 


The second of these is remarkable in that it contains no components of a. 
Instead, it provides a condition on f. In mechanical terms, it states that 
there can be no equilibrium unless the horizontal forces on the two masses 
are equal and opposite. Of course one could have observed this directly 
from the layout of the truss. In modern, three--dimensional structures with 
thousands of members meant to shelter or convey humans one should not 
however be satisfied with the ‘visual’ integrity of the structure. In particular, 
one desires a detailed description of all loads that can, and, especially, all 
loads that can not, be equilibrated by the proposed truss. In algebraic terms, 
given a matrix S, one desires a characterization of 


1. all those f for which Sx = f possesses a solution 
2. all those f for which Sa = f does not possess a solution 


We will eventually provide such a characterization in our later discussion of 
the column space of a matrix. 


Supposing now that f; + f3 = 0 we note that although the system above is 
consistent it still fails to uniquely determine the four components of z. In 
particular, it specifies only the difference between x; and x3. As a result 
both 


fi 0 
ko 
j fe 
12 ky 
2=|™ |anda = A 
0 kee 
fa fa 
kg k3 


satisfy Sa = f. In fact, one may add to either an arbitrary multiple of 


Equation: 


N 
I 
oF oO RF 


and still have a solution of Sa = f. Searching for the source of this lack of 
uniqueness we observe some redundancies in the columns of S. In 
particular, the third is simply the opposite of the first. As S' is simply 
A'’KAwe recognize that the original fault lies with A, where again, the 
first and third columns are opposites. These redundancies are encoded in z 
in the sense that 


Az=0 


Interpreting this in mechanical terms, we view z as a displacement and Az 
as the resulting elongation. In 


Az=0 


we see a nonzero displacement producing zero elongation. One says in this 
case that the truss deforms without doing any work and speaks of z as an 
unstable mode. Again, this mode could have been observed by a simple 
glance at [link]. Such is not the case for more complex structures and so the 
engineer seeks a systematic means by which all unstable modes may be 
identified. We shall see later that all these modes are captured by the null 
space of A. 


From 
S2=0 


one easily deduces that S' is singular. More precisely, if S~+ were to exist 
then S~!S'z would equal S$ 10, i.e, z =O, contrary to [link]. As a result, 
Matlab will fail to solve Sa = f even when f is a force that the truss can 


equilibrate. One way out is to use the pseudo-inverse, as we shall see in the 
General Planar Truss module. 


Matrix Methods for Mechanical Systems: The General Planar Truss 


Let us now consider something that resembles the mechanical prospection 
problem introduced in the introduction to matrix methods for mechanical 
systems. In the figure below we offer a crude mechanical model of a planar 
tissue, Say, e.g., an excised sample of the wall of a vein. 

A crude tissue model 
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Elastic fibers, numbered 1 through 20, meet at nodes, numbered 1 through 9. 
We limit our observation to the motion of the nodes by denoting the horizontal 
and vertical displacements of node 7 by x2;_1 (horizontal) and x9; (vertical), 
respectively. Retaining the convention that down and right are positive we note 
that the elongation of fiber 1 is 


€, = 22 — XB 
while that of fiber 3 is 
€3 = @3— 1. 


As fibers 2 and 4 are neither vertical nor horizontal their elongations, in terms 
of nodal displacements, are not so easy to read off. This is more a nuisance 
than an obstacle however, for noting our discussion of elongation in the small 
planar truss module, the elongation is approximately just the stretch along its 
undeformed axis. With respect to fiber 2, as it makes the angle —7 with 


respect to the positive horizontal axis, we find 


( *) in ( *) Lg — X1 + L2— X10 
€) = @gcos| —— } — zj9sin{ —— } = ————__—_. 
4 4 V/22 


Similarly, as fiber 4 makes the angle — 3m with respect to the positive 
horizontal axis, its elongation is 


e Zz cos( T) x sin( 7) ik a a Par) 
4= 27 —7> |] — & 4 
4 4 /92 


These are both direct applications of the general formula 
Equation: 


€; = Ln—-1 C08(8;) — Lan sin(O;) 


for fiber 7, as depicted in [link] below, connecting node m to node n and 
making the angle 6; with the positive horizontal axis when node m is assumed 
to lie at the point (0,0). The reader should check that our expressions for e; and 
e3 indeed conform to this general formula and that e2 and e, agree with ones 
intuition. For example, visual inspection of the specimen suggests that fiber 2 
can not be supposed to stretch (i.e., have positive eg) unless zg > x, and/or 

£2 > X19. Does this jive with [link]? 


original 


deformed 


Elongation of a generic bar, see [link]. 


Applying [link] to each of the remaining fibers we arrive ate = Aw where A 
is 20-by-18, one row for each fiber, and one column for each degree of 
freedom. For systems of such size with such a well defined structure one 
naturally hopes to automate the construction. We have done just that in the 


accompanying M-file and diary. The M-file begins with a matrix of raw data 
that anyone with a protractor could have keyed in directly from [link]: 


data = [ % 
one row of data for each fiber, the 

1 4 -pi/2 % first two 
columns are starting and ending 

1 5 -pi/4 % node 
numbers, respectively, while the third is the 

1 2 0 % angle the 
fiber makes with the positive horizontal axis 

2 4 -3*pi/4 

...and so on... ] 


This data is precisely what [link] requires in order to know which columns of 
A receive the proper cos or sin. The final A matrix is displayed in the diary. 


The next two steps are now familiar. If kK denotes the diagonal matrix of fiber 
stiffnesses and f denotes the vector of nodal forces then 


y= Ke and A'y=f 


and so one must solve Sa = f where S = A’? KA. In this case there is an 
entire three--dimensional class of z for which Az = 0 and therefore Sz = 0. 
The three indicates that there are three independent unstable modes of the 
specimen, e.g., two translations and a rotation. As a result S is singular and x 
= S\f in MATLAB will get us nowhere. The way out is to recognize that S 
has 18 — 3 = 15 stable modes and that if we restrict S to 'act' only in these 
directions then it “should' be invertible. We will begin to make these notions 
precise in discussions on the Fundamental Theorem of Linear Algebra. For 
now let us note that every matrix possesses such a pseudo-inverse and that it 
may be computed in MATLAB via the pinv command. Supposing the fiber 
stiffnesses to each be one and the edge traction to be of the form 


f=(-1 10111 =100010 =) +1 0 -1 1 =1), 


we alrive at x via X=pinv(S)*f and offer below its graphical 
representation. 


Before-After Plot 


Before and after shots of 
the truss in [link]. The 
solid (dashed) circles 

correspond to the nodal 
positions before (after) 
the application of the 

traction force, f. 


Exercises for Matrix Methods for Mechanical Systems 
Exercise: 


Problem: With regard to the uniaxial truss figure, 


¢ (i) Derive the A and K matrices resulting from the removal of the 
fourth spring, 

¢ (ii) Compute the inverse, by hand via Gauss-Jordan, of the 
resulting A? KA with k, =~ k, =k, =k 

e (iii) Use the result of (ii) to find the displacement corresponding 
to the load f = (00F)’. 


Exercise: 


Problem: 


connected by 42 fibers. Introduce one stiff (say & = 100) fiber and 
show how to detect it by 'properly' choosing f. Submit your well- 
documented M-file as well as the plots, similar to the before-after plot 
in the general planar module, from which you conclude the presence of 
a stiff fiber. 


A copy of the before-after 
figure from the general 
planar module. 


Column Space 


The Column Space 


We begin with the simple geometric interpretation of matrix-vector 
multiplication. Namely, the multiplication of the n-by-1 vector x by the m- 
by-n matrix A produces a linear combination of the columns of A. More 
precisely, if a; denotes the 7th column of A, then 

Equation: 


Ar = (a, a... Gn) 


Y1Q1, + %9A9 +... + 27,ay, 


The picture that I wish to place in your mind's eye is that Az lies in the 
subspace spanned by the columns of A. This subspace occurs so frequently 
that we find it useful to distinguish it with a definition. 


Column Space 
The column space of the m-by-n matrix S is simply the span of the its 
columns, i.e. Ra(S) = {Sx|x € R"}. This is a subspace of SR. The 
notation Ra stands for range in this context. 


Example 


Let us examine the matrix: 
Equation: 


The column space of this matrix is: 


Equation: 
0 1 0 0 
Ra(A)= 2, —-1l +2. 0 +23 1 +a, 0 2xeER?* 
0 0 0 1 


As the third column is simply a multiple of the first, we may write: 
Equation: 


0 | 0 
Ra(A)= 2 1 +22 0 +23 0 «ER? 
0 0 a) 


As the three remaining columns are linearly independent we may go no 
further. In this case, Ra (A) comprises all of R°. 


Method for Finding a Basis 


To determine the basis for Ra (A) (where A is an arbitrary matrix) we must 
find a way to discard its dependent columns. In the example above, it was 
easy to see that columns 1 and 3 were colinear. We seek, of course, a more 
systematic means of uncovering these, and perhaps other less obvious, 
dependencies. Such dependencies are more easily discerned from the row 
reduced form. In the reduction of the above problem, we come very easily 
to the matrix 

Equation: 


0 
Ared = 0 1 
0 0 


oo 
- Oo © 


Once we have done this, we can recognize that the pivot column are the 
linearly independent columns of A;-g. One now asks how this might help us 
distinguish the independent columns of A. For, although the rows of Aye 
are linear combinations of the rows of A, no such thing is true with respect 
to the columns. The answer is: pay attention to the indices of the pivot 
columns. In our example, columns {1, 2, 4} are the pivot columns of Ayeg 
and hence the first, second, and fourth columns of A, i.e., 

Equation: 


comprise a basis for Ra (A). In general: 


A Basis for the Column Space 
Suppose A is m-by-n. If columns {c,| 7 = 1, ...,r} are the pivot 
columns of Ayeq then columns {c;| j = 1,...,r} of A constitute a basis 
for Ra (A). 


Null Space 


Null Space 
Null Space 


The null space of an m-by-n matrix A is the collection of those vectors 
in R” that A maps to the zero vector in IR”. More precisely, 


N (A) = {x € R”| Ax = 0} 


Null Space Example 


As an example, we examine the matrix A: 
Equation: 


It is fairly easy to see that the null space of this matrix is: 
Equation: 


tEeR 


~ 
a 
— 
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Cs 
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This is a line in R*. 


The null space answers the question of uniqueness of solutions to Sx = f. 
For, if Sx = f and Sy = f then S(a — y) = Sw —-Sy=f—f=0 
and so (a — y) € V(S). Hence, a solution to Sa = f will be unique if, 
and only if, W(S) = {0}. 


Method for Finding the Basis 


Let us now exhibit a basis for the null space of an arbitrary matrix A. We 
note that to solve Aw = 0 is to solve A,.qa = 0. With respect to the latter, 
we suppose that 

Equation: 


{c;| j= {1, tesy rt} 


are the indices of the pivot columns and that 
Equation: 


{e;|j = {r+1,...,n}} 


are the indices of the nonpivot columns. We accordingly define the r pivot 
variables 
Equation: 


eee j= {1, a rh} 


and the n — r free variables 
Equation: 


ee i {r+1,...,n}} 


One solves A,-ga = 0 by expressing each of the pivot variables in terms of 
the nonpivot, or free, variables. In the example above, 21, %2, and x4 are 
pivot while 23 is free. Solving for the pivot in terms of the free, we find 

x4 = 0, x3 = 1, X2 = DO, or, written as a vector, 

Equation: 


where 3 is free. As x3 ranges over all real numbers the x above traces out a 
line in R*. This line is precisely the null space of A. Abstracting these 
calculations we arrive at: 


A Basis for the Null Space 
Suppose that A is m-by-n with pivot indices {c,| 7 = {1,...,r}} and 
free indices {c,| 7 = {r+1,...,n}}. A basis for (A) may be 
constructed of n — r vectors {2', z?,...,z”-"} where z*, and only z*, 
possesses a nonzero in its c,4;, Component. 


A MATLAB Observation 


As usual, MATLAB has a way to make our lives simpler. If you have 
defined a matrix A and want to find a basis for its null space, simply call the 
function nul1(A). One small note about this function: if one adds an extra 
flag, 'r',asinnull(A, 'r' ), then the basis is displayed "rationally" 
as opposed to purely mathematically. The MATLAB help pages define the 
difference between the two modes as the rational mode being useful 
pedagogically and the mathematical mode of more value (gasp!) 
mathematically. 


Final thoughts on null spaces 


There is a great deal more to finding null spaces; enough, in fact, to warrant 
another module. One important aspect and use of null spaces is their ability 
to inform us about the uniqueness of solutions. If we use the column space 
to determine the existence of a solution a to the equation Aw = b. Once we 
know that a solution exists it is a perfectly reasonable question to want to 
know whether or not this solution is the only solution to this problem. The 


hard and fast rule is that a solution x is unique if and only if the null space 
of A is empty. One way to think about this is to consider that if Aw = 0 
does not have a unique solution then, by linearity, neither does Aw = b. 
Conversely, if (Az = 0) A (240) A (Ay=b) then A (z+ y) = bas 
well. 


The Null Space and Column Space: An Example 


Preliminary Information 


Let us compute bases for the null and column spaces of the adjacency 
matrix associated with the ladder below. 
An unstable ladder? 


The ladder has 8 bars and 4 nodes, so 8 degrees of freedom. Denoting the 
horizontal and vertical displacements of node 7 by x2;_1 and 72;, 
respectively, we arrive at the A matrix 


1 0 0 00 0 0 
=. 0 1 /0 © ©: 0 0 
0 0 -1 0 00 0 0 
0 =1 0 0 0 1 0 0 
ee Gy ek “Oe OE 
0 0 0 0 10 0 0 
6: <0: oO 0: Oe = 0 
0 0 0 0 0 0-1 0 


Finding a Basis for the Column Space 


To determine a basis for #(A) we must find a way to discard its dependent 
columns. A moment's reflection reveals that columns 2 and 6 are colinear, 
as are columns 4 and 8. We seek, of course, a more systematic means of 


uncovering these and perhaps other less obvious dependencies. Such 
dependencies are more easily discerned from the row reduced form 


10000 00 0 
(61000 -1.0 6 
00100000 
00010 0 0-1 

Be NAS Foy a. GHG 
0000001 0 
0000000 0 
0000000 0 


Recall that rref performs the elementary row operations necessary to 
eliminate all nonzeros below the diagonal. For those who can't stand to miss 
any of the action I recommend rrefmovie. 


Each nonzero row of A,.q is called a pivot row. The first nonzero in each 
row of A,.g is called a pivot. Each column that contains a pivot is called a 
pivot column. On account of the staircase nature of A,,g we find that there 
are aS many pivot columns as there are pivot rows. In our example there are 
six of each and, again on account of the staircase nature, the pivot columns 
are the linearly independent columns of A,.4g. One now asks how this might 
help us distinguish the independent columns of A. For, although the rows of 
Ayeq are linear combinations of the rows of A, no such thing is true with 
respect to the columns. In our example, columns {1, 2,3, 4, 5,7} are the 
pivot columns. In general: 

Proposition 


Suppose A is m-by-n. If columns {c,| 7 = 1, ...,r} are the pivot columns of 
Area. then columns {c,| 7 = 1,...,r} of A constitute a basis for (A). 


Note that the pivot columns of A;eg. are, by construction, linearly 
independent. Suppose, however, that columns {c,| 7 = 1,...,r} of A are 


linearly dependent. In this case there exists a nonzero x € IR” for which 
Az = 0 and 
Equation: 


Vek C46|j = lant} t @e—0) 


Now Aa = O necessarily implies that A,.ga = 0, contrary to the fact that 
columns {c;| 7 = 1, ...,r} are the pivot columns of Area. 


We now show that the span of columns {c,| 7 = 1,...,r} of A indeed 
coincides with &(A). This is obvious if r = n, ie., if all of the columns 
are linearly independent. If r < n, there exists ag ¢ {c,|j = 1,...,r}. 
Looking back at A,,g. we note that its gth column is a linear combination of 
the pivot columns with indices not exceeding g. Hence, there exists an x 
satisfying [link] and A;.qa = O, and x, = 1. This z then necessarily 
satisfies Aw = O. This states that the gth column of A is a linear 
combination of columns {c,| 7 = 1,...,r} of A. 


Finding a Basis for the Null Space 


Let us now exhibit a basis for (A). We exploit the already mentioned 
fact that. VW (A) = W(Ayjea). Regarding the latter, we partition the 
elements of x into so called pivot variables, 


1 ees ee F ra, 
and free variables 
{zz|k ¢ {c,|7 = 1,...,r}} 


There are evidently n — r free variables. For convenience, let us denote 
these in the future by 


eee g=f-—-1, il 


One solves A,-ga = O, by expressing each of the pivot variables in terms 
of the nonpivot, or free, variables. In the example above, x1, £9, ©3, £4, Zs, 
and 2x7 are pivot while xg and xg are free. Solving for the pivot in terms of 
the free we find 


L7 = 0 
5 = 0 
“LA = LB 
L3 = 0 
“2 = 6 
Li = 0 
or, written as a vector, 
Equation: 
0 0 
1 0) 
0 0 
0 1 
= 7% 0 + £8 0 
1 0 
0 0 
0 1 


where 2g and xg are free. As xg and xg range over all real numbers, the a 
above traces out a plane in R®. This plane is precisely the null space of A 
and [link] describes a generic element as the linear combination of two 
basis vectors. Compare this to what MATLAB returns when faced with 
null(A, 'r' ). Abstracting these calculations we arrive at 

Proposition 


Suppose that A is m-by-n with pivot indices {c,| 7 = 1,...,r} and free 
indices {c,|j =r+1,...,n}. A basis for W(A). may be constructed of 

n — r vectors ee BP sts ety where z*, and only zk possesses a nonzero 
in its c,4, Component. 


The Physical Meaning of Our Calculations 


Let us not end on an abstract note however. We ask what #(A) and (A) 
tell us about the ladder. Regarding #(A) the answer will come in the next 
chapter. The null space calculation however has revealed two independent 
motions against which the ladder does no work! Do you see that the two 
vectors in [link] encode rigid vertical motions of bars 4 and 5 respectively? 
As each of these lies in the null space of A, the associated elongation is 
zero. Can you square this with the ladder as pictured in [link]? I hope not, 
for vertical motion of bar 4 must 'stretch' bars 1, 2, 6, and 7. How does one 
resolve this (apparent) contradiction? 


Left Null Space 


Left Null Space 


If one understands the concept of a null space, the left null space is 
extremely easy to understand. 


Left Null Space 
The Left Null Space of a matrix is the null space of its transpose, i.e., 


NW (AT) ={yeR™ ATy=0} 


The word "left" in this context stems from the fact that A? y = 0 is 
equivalent to y?.A = 0 where y "acts" on A from the left. 


Example 


As Ayeq was the key to identifying the null space of A, we shall see that 
AZ , is the key to the null space of A’. If 


Equation: 
1 1 
As A. 2 
1 3 
then 
Equation: 
AP 111 
t23 
and so 


Equation: 


111 
A= 
red é 1 ,) 


We solve ACs = 0 by recognizing that y; and y2 are pivot variables while 
y3 is free. Solving Aly = 0 for the pivot in terms of the free we find 


Sp (2y3) and y; = y3 hence 
Equation: 


1 
N(AT)= ys -2 yzsER 
1 


Finding a Basis for the Left Null Space 


The procedure is no different than that used to compute the null space of A 
itself. In fact 


A Basis for the Left Null Space 
Suppose that A” is n-by-m with pivot indices {c,;| 7 = {1,...,r}} 
and free indices {c;| 7 = {r+1,..., m}}. A basis for. (A*) may 
be constructed of m — r vectors ea en ge) where z*, and 
only 2k possesses a nonzero in its c,4;% Component. 


Row Space 


The Row Space 


As the columns of A? are simply the rows of A we call Ra (A*) the row 
space of A? . More precisely 


Row Space 
The row space of the m-by-n matrix A is simply the span of its rows, 
Le. 
Ra (A‘) = {A’y ye R™} 


This is a subspace of IR”. 


Example 


Let us examine the matrix: 


Equation: 
0 10 0 
A= -101 0 
0 00 1 
The row space of this matrix is: 
Equation: 
0 =] 0 
1 0 0 
Ra(A*)= y + yp + y3 yeR 
0 | 0 
0 0 1 


As these three rows are linearly independent we may go no further. We 
"recognize" then #a(A7) as a three dimensional subspace of R*. 


Method for Finding the Basis of the Row Space 


Regarding a basis for Ba (A*) we recall that the rows of Ayeq, the row 


reduced form of the matrix A, are merely linear combinations of the rows 
of A and hence 
Equation: 


Ra (A*) = Sa Ayea) 


This leads immediately to: 


A Basis for the Row Space 
Suppose A is m-by-n. The pivot rows of A;eq constitute a basis for 
Ra (A*) ; 


With respect to our example, 
Equation: 


— 
_ 
Foo;o 


comprises a basis for Za (A7). 


Exercises: Columns and Null Spaces 
Exercises 


1. I encourage you to use rref and nul for the following. 


o (i) Add a diagonal crossbar between nodes 3 and 2 in the unstable 
ladder figure and compute bases for the column and null spaces of 
the new adjacency matrix. As this crossbar fails to stabilize the 
ladder, we shall add one more bar. 

o (ii) To the 9 bar ladder of (i) add a diagonal cross bar between 
nodes 1 and the left end of bar 6. Compute bases for the column 
and null spaces of the new adjacency matrix. 


2. We wish to show that N(.A) = N(A™A) regardless of A. 


o (i) We first take a concrete example. Report the findings of nul1 
when applied to A and A? A for the A matrix associated with the 
unstable ladder figure. 

o (ii) Show that N(A) C N(A7A), ie. that if Aw = 0 then 
A’ Ax = 0. 

o (iii) Show that N(A7A) C N(A), ie., that if A’ Aa = 0 then 
Az = 0. (Hint: if A7Ax = 0 then x7 A? Az = 0.) 


3. Suppose that A is m-by-n and that V(A) = R”. Argue that A must be 
the zero matrix. 


Vector Space 


Introduction 


You have long taken for granted the fact that the set of real numbers, R, is 
closed under addition and multiplication, that each number has a unique 
additive inverse, and that the commutative, associative, and distributive 
laws were right as rain. The set, C, of complex numbers also enjoys each of 
these properties, as do the sets IR” and C” of columns of n real and complex 
numbers, respectively. 


To be more precise, we write a and y in R” as 
= T 
SS pire<, ea) 


T 
y = (y1Yy2- - -Yn) 


and define their vector sum as the elementwise sum 
Equation: 


Lit yi 


2+ Y2 
e+y= 


Ln + Yn 


and similarly, the product of a complex scalar, z € C with @ as: 
Equation: 


Vector Space 


These notions lead naturally to the concept of vector space. A set V is said 
to be a vector space if 


le+y=y+2 foreach z and yinV 

2e+y+z=x2+y+ 2 foreach ag, y, and zinV 

3. There is a unique "zero vector" such that 2 + O = a for each x in V 
4. For each x in V there is a unique vector —2 such that 2 + —a2 = O. 
5. Lae a 

6. (cyc2)@ = cy (cpx) for each w in V and cy and cp in C. 

7.c(a@ + y) = ca + cy for each w and y in V andcinC. 

8. (cy + c2)@ = cyx@ + Cow for each w in V and c; and cz in C. 


Subspaces 


Subspace 


A subspace is a subset of a vector space that is itself a vector space. The 
simplest example is a line through the origin in the plane. For the line is 
definitely a subset and if we add any two vectors on the line we remain on 
the line and if we multiply any vector on the line by a scalar we remain on 
the line. The same could be said for a line or plane through the origin in 3 
space. As we shall be travelling in spaces with many many dimensions it 
pays to have a general definition. 


A subset S' of a vector space V is a subspace of V when 
if x and y belong to S then so does x + y. 
if x belongs to S and ¢ is real then tz belong to S. 


As these are oftentimes unwieldy objects it pays to look for a handful of 
vectors from which the entire subset may be generated. For example, the set 
of x for which x; + 22 + #3 + 24 = 0 constitutes a subspace of R*. Can 
you 'see' this set? Do you 'see' that 


and 


and 


not only belong to a set but in fact generate all possible elements? More 
precisely, we say that these vectors span the subspace of all possible 
solutions. 


Span 
A finite collection {s1, 52,...,8,} of vectors in the subspace S is said 
to span S if each element of S can be written as a linear combination 
of these vectors. That is, if foreach s € S there exist n reals 
{21, £2,..-,£n} such that s = 118] + £989 +... + EnSn. 


When attempting to generate a subspace as the span of a handful of vectors 
it is natural to ask what is the fewest number possible. The notion of linear 


independence helps us clarify this issue. 


Linear Independence 


A finite collection {s1, 2,..., 8} of vectors is said to be linearly 
independent when the only reals, {21, £2,..., 2} for which 
%+%+...+ 2, =Oarexr; = 2%. =...= 2, = DO. In other words, 


when the null space of the matrix whose columns are {81, 52,..., Sn} 
contains only the zero vector. 


Combining these definitions, we arrive at the precise notion of a ‘generating 
set.’ 


Basis 


Any linearly independent spanning set of a subspace S is called a basis 
of S. 


Though a subspace may have many bases they all have one thing in 
common: 


Dimension 


The dimension of a subspace is the number of elements in its basis. 


Row Reduced Form 


Row Reduction 


A central goal of science and engineering is to reduce the complexity of a 
model without sacrificing its integrity. Applied to matrices, this goal 
suggests that we attempt to eliminate nonzero elements and so 'uncouple' 
the rows. In order to retain its integrity the elimination must obey two 
simple rules. 


Elementary Row Operations 
You may swap any two rows. 
You may add to a row a constant multiple of another row. 


With these two elementary operations one can systematically eliminate all 
nonzeros below the diagonal. For example, given 
Equation: 


an) 

no Oo oO F 
wonorF © 
ae © CO 


it seems wise to swap the first and fourth rows and so arrive at 
Equation: 


1234 
0 1 0 0 
—1 0 1 0 
0 00 1 


adding the first row to the third now produces 
Equation: 


123 4 
0 10 0 
02 4 4 
00 0 1 


subtracting twice the second row from the third yields 


Equation: 
12 3 4 
010 0 
00 4 4 
0001 


a matrix with zeros below its diagonal. This procedure is not restricted to 
Square matrices. For example, given 


Equation: 
111i 
244 2 
3.5 5 3 


we Start at the bottom left then move up and right. Namely, we subtract 3 
times the first row from the third and arrive at 


Equation: 
1ii1ii1 
244 2 
02 2 0 


and then subtract twice the first row from the second, 
Equation: 


oO 
| NO 
| NO 
oO 


and finally subtract the second row from the third, 


Equation: 
111i 
02 2 0 
00 0 0 


It helps to label the before and after matrices. 


The Row Reduced Form 
Given the matrix A we apply elementary row operations until each 
nonzero below the diagonal is eliminated. We refer to the resulting 
matrix aS Ayeg. 


Uniqueness and Pivots 


As there is a certain amount of flexibility in how one carries out the 
reduction it must be admitted that the reduced form is not unique. That is, 
two people may begin with the same matrix yet arrive at different reduced 
forms. The differences however are minor, for both will have the same 
number of nonzero rows and the nonzeros along the diagonal will follow 
the same pattern. We capture this pattern with the following suite of 
definitions, 


Pivot Row 
Each nonzero row of Ayeg is called a pivot row. 


Pivot 
The first nonzero term in each row of Ayeq is called a pivot. 


Pivot Column 


Each column of A;eq that contains a pivot is called a pivot column. 


Rank 
The number of pivots in a matrix is called the rank of that matrix. 


Regarding our example matrices, the first has rank 4 and the second has 
rank 2. 


Row Reduction in MATLAB 


MATLAB's rref command goes full-tilt and attempts to eliminate ALL 
off diagonal terms and to leave nothing but ones on the diagonal. I 
recommend you try it on our two examples. You can watch its individual 
decisions by using rrefmovie instead. 


Least Squares 


Introduction 


We learned in the previous chapter that Ax = 6 need not possess a solution 
when the number of rows of A exceeds its rank, i.e., 7 < m. As this 
situation arises quite often in practice, typically in the guise of 'more 
equations than unknowns,' we establish a rationale for the absurdity Ax = b 


The Normal Equations 


The goal is to choose x such that Az is as close as possible to b. Measuring 
closeness in terms of the sum of the squares of the components we arrive at 
the 'least squares’ problem of minimizing 
Equation: 

res 


(|| Az —b ||)” = (Aw — 6)" (Aw — 8) 


over all x € R. The path to the solution is illuminated by the Fundamental 
Theorem. More precisely, we write 
Vor, by, be € R(A) A bn € N(A’*) : (b = bey + by). On noting that (i) 
Vbp, x € R”: ((Az — br) € R(A)) and (ii)R(A) L N(A*) we arrive at 
the Pythagorean Theorem. 
Equation: 

Pythagorean Theorem 


norm?(Az —b) = (|| Aw —bg+by |\)’ 
(|| Az — bp |)” + (|| by Il)” 


| 


It is now clear from the Pythagorean Theorem that the best z is the one that 
satisfies 
Equation: 


Az = br 


As br € R(A) this equation indeed possesses a solution. We have yet 
however to specify how one computes bz given b. Although an explicit 
expression for bp, the so called orthogonal projection of b onto R(A), in 
terms of A and 0 is within our grasp we shall, strictly speaking, not require 
it. To see this, let us note that if x satisfies the above equation then 
Equation: 


Az—b = Ax—brt+ bn 


—bn 


As by is no more easily computed than br you may claim that we are just 
going in circles. The 'practical' information in the above equation however 
is that (Ax — b) € A?,ie., A? (Ax — b) = 0,ie., 

Equation: 


A’ Ax = A™b 


As AlbE R(A*) regardless of b this system, often referred to as the 
normal equations, indeed has a solution. This solution is unique so long as 
the columns of A”’A are linearly independent, i.e., so long as 

N(A‘A) = {0}. Recalling Chapter 2, Exercise 2, we note that this is 


equivalent to N(A) = {0}. We summarize our findings in 


The set of « € br for which the misfit (|| Aa — 6 ||)? is smallest is 
composed of those x for which A? Ax = A‘b. There is always at least one 
such x. There is exactly one such x if N(A) = {0}. 


As a concrete example, suppose with reference to [link] that A = 


oo 
oe 


andb=)1 


The decomposition of 6. 


As b # R(A) there is no x such that Ax = b. Indeed, 
(|| Ax — b ||)? = (#1 + w+ —1)* + (x, — 1)? +1 > 1, with the 


minimum uniquely attained at x = i}? in agreement with the unique 


1 1 1 
solution of the above equation, for ATA = ( ;) and A’b = () _We 


now recognize, a posteriori, that bk = Ax = | 1 | is the orthogonal 


0 
projection of b onto the column space of A. 


Applying Least Squares to the Biaxial Test Problem 


We shall formulate the identification of the 20 fiber stiffnesses in this 
previous figure, as a least squares problem. We envision loading, f, the 9 
nodes and measuring the associated 18 displacements, x. From knowledge 
of x and f we wish to infer the components of kK = diag (k) where k is the 
vector of unknown fiber stiffnesses. The first step is to recognize that 
Equation: 


A'K Az =f 


may be written as 
Equation: 


VB, B = A’ diag (Az) : (Bk = f) 


Though conceptually simple this is not of great use in practice, for B is 18- 
by-20 and hence the above equation possesses many solutions. The way out 
is to compute k as the result of more than one experiment. We shall see that, 
for our small sample, 2 experiments will suffice. To be precise, we suppose 
that x1 is the displacement produced by loading f! while x? is the 
displacement produced by loading f?. We then piggyback the associated 


(Al diag (Az) P\ 
pieces in B = —_ x, } and f= , } This B is 36-by-20 and 
A’ diag (Ax ) f 


so the system Bk = f is overdetermined and hence ripe for least squares. 


We proceed then to assemble B and f. We suppose f! and f? to correspond 
to horizontal and vertical stretching 
Equation: 


fi=(-100010-100010-10001 0) 
Equation: 


fP=0101010000000 -1 0 -1 0 -1)* 


respectively. For the purpose of our example we suppose that each k; = 1 
except kg = 5. We assemble A’ K A as in Chapter 2 and solve 
Equation: 


A’ K Az! = f 


with the help of the pseudoinverse. In order to impart some ‘reality’ to this 
problem we taint each x with 10 percent noise prior to constructing B. 
Regarding 

Equation: 


B' Bk = B'f 
we note that Matlab solves this system when presented with kK=B\f when 


B is rectangular. We have plotted the results of this procedure in the [Link]. 
The stiff fiber is readily identified. 


fiber stiffness 
to h& 


ho 


15 20 25 


10 
fiber number 


Results of a successful biaxial test. 


Projections 


From an algebraic point of view [link])is an elegant reformulation of the 
least squares problem. Though easy to remember it unfortunately obscures 
the geometric content, suggested by the word ‘projection,’ of [link]. As 
projections arise frequently in many applications we pause here to develop 
them more carefully. With respect to the normal equations we note that if 
N(A) = {0} then 

Equation: 


z= (ATA) ‘AT 


and so the orthogonal projection of b onto R(A) is: 


Equation: 
br — Ax 
A(ATA) *ATD 
Defining 
Equation: 


P= A(ATA) ‘AT 


[link] takes the form bp = Pb. Commensurate with our notion of what a 
‘projection’ should be we expect that P map vectors not in R(A) onto R(A) 
while leaving vectors already in IR(A) unscathed. More succinctly, we 
expect that Pbr = bp, i.e., Phy = Por. As the latter should hold for all 

b € R™ we expect that 

Equation: 


P?=P 


With respect to [link] we find that indeed 


Equation: 
P? = A(ATA)“ATA(ATA) “AT 
A(ATA) "AT 
= Pp 


We also note that the P in [link] is symmetric. We dignify these properties 
through 


orthogonal projection 
A matrix P that satisfies P? = P is called a projection. A symmetric 
projection is called an orthogonal projection. 


We have taken some pains to motivate the use of the word ‘projection.’ You 
may be wondering however what symmetry has to do with orthogonality. 
We explain this in terms of the tautology 

Equation: 


b= Pb— Ib 


Now, if P is a projection then so too is f — P. Moreover, if P is symmetric 
then the dot product of b's two constituents is 


Equation: 

(Pb) (I-— P)b = b'PT(I-P)b 
= b'(P-P’*)b 
= b'0b 

0 


i.e., Pb is orthogonal to (I — P)b. As examples of a nonorthogonal 


1 0 O 
projections we offer P = — O OJ] andJ — P. Finally, let us note 
av: et 
ye ae 


that the central formula, P = A(A™A) ee A’, is even a bit more general 


than advertised. It has been billed as the orthogonal projection onto the 
column space of A. The need often arises however for the orthogonal 
projection onto some arbitrary subspace M/. The key to using the old P is 
simply to realize that every subspace is the column space of some matrix. 
More precisely, if 

Equation: 


aCe nee 


is a basis for M then clearly if these x; are placed into the columns of a 
matrix called A then R(A) = M. For example, if M is the line through 
(1 1)" then 

Equation: 


is orthogonal projection onto M. 


Exercises 


1. Gilbert Strang was stretched on a rack to lengths £ = 6, 7, and 8 feet 
under applied forces of f = 1, 2, and 4 tons. Assuming Hooke's law 
£ — L = cf, find his compliance, c, and original height, L, by least 
squares. 

2. With regard to the example of § 3 note that, due to the the random 
generation of the noise that taints the displacements, one gets a 
different 'answer' every time the code is invoked. 


1. Write a loop that invokes the code a statistically significant 
number of times and submit bar plots of the average fiber 
stiffness and its standard deviation for each fiber, along with the 
associated M--file. 


2. Experiment with various noise levels with the goal of determining 


the level above which it becomes difficult to discern the stiff 
fiber. Carefully explain your findings. 


3. Find the matrix that projects R® onto the line spanned by (1 0 ioe 

A, Find the matrix that projects IR® onto the plane spanned by (1 0 1 
and (la) 

5. If P is the projection of IR” onto a k--dimensional subspace M, what 
is the rank of P and what is R(P)? 


Nerve Fibers and the Dynamic Strang Quartet 


Introduction 
Up to this point we have largely been concerned with 
1. Deriving linear systems of algebraic equations (from considerations of static 
equilibrium) and 


2. The solution of such systems via Gaussian elimination. 


In this module we hope to begin to persuade the reader that our tools extend in a natural 
fashion to the class of dynamic processes. More precisely, we shall argue that 


1. Matrix Algebra plays a central role in the derivation of mathematical models of 
dynamical systems and that, 


2. With the aid of the Laplace transform in an analytical setting or the Backward Euler 
method in the numerical setting, Gaussian elimination indeed produces the solution. 


Nerve Fibers and the Dynamic Strang Quartet 


Gathering Information 


A nerve fiber's natural electrical stimulus is not direct current but rather a short burst of 
current, the so-called nervous impulse. In such a dynamic environment the cell's membrane 
behaves not only like a leaky conductor but also like a charge separator, or capacitor. 

An RC model of a nerve fiber 


where pF’ denotes micro-Farad. Recalling our variable conventions, the capacitance of a 
single compartment is 


l 
Cy, = 27a—c 
N 


and runs parallel to each R,,, see [link]. This figure also differs from the simpler circuit from 
the introductory electrical modeling module in that it possesses two edges to the left of the 
stimuli. These edges serve to mimic that portion of the stimulus current that is shunted by the 
cell body. If A.p denotes the surface area of the cell body, then it has 


capacitance of cell body 
Cob = Acbe 


resistance of cell body 
Re = AcbPm: 


Updating the Strang Quartet 


We ask now how the static Strang Quartet of the introductory electrical module should be 
augmented. 


Updating (S1') 
Regarding (S1') we proceed as before. The voltage drops are 
ey = ZX 
eg = ®t Ey. 
€3 = %— r 
€4 = X22 
€5 = @2— Ee 
€6 = %2—- 2X3 
C7 = 23 
€3g = &3— pan 


and so 


0 —1 0 0 
1 —1 0 0 
0 —1 1 0 
0 0 -1 O 
e=b-—Azx where b= (-E,,) i and A= . <f. @ 
0 0 -1 1 
0 0 0 -l1 
1 0 Q -l 


Updating (S2) 
To update (S2) we must now augment Ohm's law with 


Voltage-current law obeyed by a capacitor 
The current through a capacitor is proportional to the time rate of change of the 
potential across it. 


This yields, (denoting derivative by '), 


Y= Copei 
€2 
= 
u iis 
¥3 R, 
ys = Cea’ 
€5 
Y= Re 
Yo = R, 
y= Cme7 
€8 
¥3 = R,, 
or, in matrix terms, 
y—=Ge+Ce’ 


where 


Q 

\| 
eco coco So oo 
a ——) 
oe 6 oo Ss 6 6 
a | Se cs 2 .S 
Sona oo So co So 
aoe oo 6 oS 
A+ STO COC CO oO 


and 


ooo oO Oo 

— ee ee ee eee ee) 
en ee ee ee) 
Sees Se 6 
oooo°eclUcUOWUlUCcCOlCUhO 
an a ee ee ee) 
oie S6-6'e°s 
ee ee ee eee) 


are the conductance and capacitance matrices. 


Updating (S3) 


As Kirchhoff's Current law is insensitive to the type of device occupying an edge, step (S3) 
proceeds exactly as before. 


o- Yi — ¥2—- y3 = 0 
¥3— Y4— Y5 — Yo = O 
¥6 — ¥7— ys = 0 
or, in matrix terms, 


i 
20 
A'y=-—f where f= 0 
0 


Step (S4): Assembling 
Step (S4) remains one of assembling, 
(ATy = —f) = (A” (Ge + Ce’) = —-f) = (A (G(b— Ax) +C (b'— Az’)) = —-f) 


becomes 
Equation: 


A'C Aa’ + A’GAa = A’Gb+ f+ A™CU’. 


This is the general form of the potential equations for an RC circuit. It presumes of the user 
knowledge of the initial value of each of the potentials, 


Equation: 
x(0) = X 
Regarding the circuit of [link], and letting G = a we find 
Cob 0 0 Gop qe Gi —G; 0 
A'CA= 0 C0 , A'TGA= —Gi  2G:+Gn  —G; 
0 0C 0 -G; -Gj+ Gm 
Gob 0 
A'Gb=Em Gm and A™Cb'= 0 
Gin 0 
and an initial (rest) potential of 
12 
#(0) =F, 1 
1 


Modes of Attack 


We shall now outline two modes of attack on such problems. The Laplace Transform is an 
analytical tool that produces exact, closed-form solutions for small tractable systems and 
therefore offers insight into how larger systems 'should' behave. The Backward-Euler 
method is a technique for solving a discretized (and therefore approximate) version of [link]. 
It is highly flexible, easy to code, and works on problems of great size. Both the Backward- 
Euler and Laplace Transform methods require, at their core, the algebraic solution of a linear 


system of equations. In deriving these methods we shall find it more convenient to proceed 
from the generic system 
Equation: 


x’ — Br+g 


With respect to our fiber problem 


Equation: 
B = (-(ATCA)')ATGA 
Cs, Cob 0 
G; —(Gi+G im) 
0 C., Ci 
and 
GoEmtio 
Co 
g=(ATCA) (ATQb+f)= 2 
EmGm 


Cm 


The Old Laplace Transform 


The Laplace Transform is typically credited with taking dynamical problems into static problems. 
Recall that the Laplace Transform of the function h is 


MATLAB is very adept at such things. For example: 


Example: 
The Laplace Transform in MATLAB 


>> syms t 

>> laplace(exp(t)) 

ans = 1/(s-1) 

>> laplace(t*(exp(-t)) 


ans = 1/(s+1)42 


The Laplace Transform of a matrix of functions is simply the matrix of Laplace transforms of the 
individual elements. 


Example: 
Laplace Transform of a matrix of functions 


module: 
x = Br+g, 


we write it as 
Equation: 


2(32) = ¥(Bx +g) 


and so must determine how -# acts on derivatives and sums. With respect to the latter it follows 
directly from the definition that 
Equation: 


L(Ba+g) (Ba) + £(g) 


= BY(x)+ L(g) 


Regarding its effect on the derivative we find, on integrating by parts, that 


oe t 
2(52) = eo) SE) ay (tye 
dt) J 


Pas f e “Ya(t) dt. 


Supposing that « and s are such that a(t)e~(**) + 0 as t > 00 we arrive at 
Equation: 


Now, upon substituting [link] and [link] into [link] we find 
sL(x) —2(0) = BY(x) + L(g), 


which is easily recognized to be a linear system for #(a), namely 
Equation: 


(sI — B) L(x) = L(g) + (0). 
The only thing that distinguishes this system from those encountered since our first brush with these 
systems is the presence of the complex variable s. This complicates the mechanical steps of Gaussian 


Elimination or the Gauss-Jordan Method but the methods indeed apply without change. Taking up the 
latter method, we write 


La) = (sI — B)* (L(g) + 2(0)). 
The matrix (sf — B) tis typically called the transfer function or resolvent, associated with B, at s. 


We turn to MATLAB for its symbolic calculation. (for more information, see the tutorial on 
MATLAB's symbolic toolbox). For example, 


Example: 


>> B = [2 -1; -1 2] 


>> R = inv(s*eye(2) -B) 
R= 
[ (s-2)/(s*s-4*s+3), -1/(s*sS-4*sS+3) ] 


[ -1/(s*s-4*s+3), (S-2)/(S*S-4*S+3) ] 


We note that (sJ — B)' is well defined except at the roots of the quadratic, s? — 4s + 3. This 


quadratic is the determinant of sf — B and is often referred to as the characteristic polynomial of B 
. Its roots are called the eigenvalues of B. 


Example: 


As a second example let us take the B matrix of the dynamic Strang quartet module with the 
parameter choices specified in fib3.m, namely 
Equation: 


—0.135 0.125 0 
[5h = >On ts 
0 0.5 —0.51 


The associated (sI — Bye: is a bit bulky (please run fib3.m) so we display here only the denominator 
of each term, i.e., 
Equation: 


s° + 1.6558” + 0.40788 + 0.0039. 


ot 
Assuming a current stimulus of the form i9(t) = vee and Em = 0 brings 


0.191 
(s+4)" 
L(g)(s)=] 9 


0 


and so [link] persists in 


soe 55 4-0.27 
0.5s + 0.26 
0.2497 


0-191 


(2) = (sl — B) L(g) = ae 
(s+ ¢) (s? + 1.655s? + 0.4078 + 0.0039) 


Now comes the rub. A simple linear solve (or inversion) has left us with the Laplace transform of z. 
The accursed 
No Free Lunch Theorem 


We shall have to do some work in order to recover # from (a). 
confronts us. We shall face it down in the Inverse Laplace module. 


The Inverse Laplace Transform 


To Come 


In The Transfer Function we shall establish that the inverse Laplace transform of a function h is 
Equation: 


il 
Or 


2mo= sf eeMn(etyit) dy 


where 7 = ¥2—1 and the real number c is chosen so that all of the singularities of h lie to the left of the line of 
integration. 


Proceeding with the Inverse Laplace Transform 


With the inverse Laplace transform one may express the solution of x’ = Ba + g, as 
Equation: 


x(t) = ¥-((s1 — B)*) (¥(g) + 2(0)) 


As an example, let us take the first component of #(a), namely 


0.19 (s? + 1.58 + 0.27) 


Ly,(s) = ; 
(s + +)*(s3 + 1.6555? + 0.4078s + 0.0039) 


We define: 


poles 
Also called singularities, these are the points s at which %;,, (s) blows up. 


These are clearly the roots of its denominator, namely 
Equation: 


o73 
—1/100, —329/400 + Se and —1/6. 


All four being negative, it suffices to take c = 0 and so the integration in [link] proceeds up the imaginary axis. We 
don't suppose the reader to have already encountered integration in the complex plane but hope that this example 
might provide the motivation necessary for a brief overview of such. Before that however we note that MATLAB 
has digested the calculus we wish to develop. Referring again to fib3.m for details we note that the ilaplace 
command produces 


/273t 
16 


a =t ~ (3294) 
a1(t) = 211.35e™ — (0.0554t° + 4.5464¢7 + 1.085t + 474.19)e* + e797 262.842 cosh 


140 rT T T T T T T 


4 n n n L 4 


100 180 200 2 300 380 400 
t (ms) 


The 3 potentials associated with the RC 
circuit model figure. 


The other potentials, see the figure above, possess similar expressions. Please note that each of the poles of #(21) 
appear as exponents in 2, and that the coefficients of the exponentials are polynomials whose degrees is 
determined by the order of the respective pole. 


The Backward-Euler Method 


Where in the Inverse Laplace Transform module we tackled the derivative 
in 
Equation: 


x’ — Br+g, 


via an integral transform we pursue in this section a much simpler strategy, 
namely, replace the derivative with a finite difference quotient. That is, one 
chooses a small dt and 'replaces' [link] with 

Equation: 


#(t) — E(t — dt) 


a = Bz(t) + 9(¢). 


The utility of [link] is that it gives a means of solving for x at the present 
time, t, from the knowledge of x in the immediate past, t — dt. For 
example, as (0) = x(0) is supposed known we write [link] as 


& - B) #(dt) = 0) + 9(dt). 


Solving this for #(dt) we return to [link] and find 


& 2 B) #(2dt) = a + g(2dt). 


and solve for %(2dt). The general step from past to present, 


Equation: 
(jdt) = (+ = B) 7 (a = alia ) 


is repeated until some desired final time, Tdt, is reached. This equation has 
been implemented in fib3.m with dt = 1 and B and g as in the dynamic 
Strang module. The resulting % ( run fib3.m yourself!) is indistinguishable 
from the plot we obtained in the Inverse Laplace module. 


Comparing the two representations, this equation and [link], we see that 
they both produce the solution to the general linear system of ordinary 
equations, see this eqn, by simply inverting a shifted copy of B. The former 
representation is hard but exact while the latter is easy but approximate. Of 
course we should expect the approximate solution, Z , to approach the exact 
solution, x, as the time step, dt , approaches zero. To see this let us return to 
[link] and assume, for now, that g = 0. In this case, one can reverse the 
above steps and arrive at the representation 

Equation: 


#(jdt) — ((Z z atB) *)'2(0) 


Now, for a fixed time t we suppose that dt = : and ask whether 


o(t) = limit ((: 2 ‘ z) : e 


This limit, at least when B is one-by-one, yields the exponential 
a(t) = e”x(0) 


clearly the correct solution to this equation. A careful explication of the 
matrix exponential and its relationship to this equation will have to wait 
until we have mastered the inverse laplace transform. 


Exercises: Matrix Methods for Dynamical Systems 


i, 


Z; 
oi 


Compute, without the aid of a machine, the Laplace transforms of et 
and te’. Show ALL of your work. 

Extract from fib3.m analytical expressions for x2 and 73. 

Use €1g to compute the eigenvalues of B as given in this equation. 
Use det to compute the characteristic polynomial of B. Use roots 
to compute the roots of this characteristic polynomial. Compare these 
to the results of €1g. How does Matlab compute the roots of a 
polynomial? (type help roots for the answer). 


. Adapt the Backward Euler portion of fib3.m so that one may specify 


an arbitrary number of compartments, as in fib1.m. Submit your 
well documented M-file along with a plot of x; and x19 versus time 
(on the same well labeled graph) for a nine compartment fiber of 
length / = 1cm. 


. Derive this equation from a previous equation by working backwards 


toward x(0). Along the way you should explain why 


=(I-a(t)B)" 


aid 
. Show, for scalar B, that (1 _ ‘B) — ePt as 7 — oo. Hint: By 


definition 


now use L'Hopital's rule to show that 7 log a — Bt. 
j 


Matrix Analysis of the Branched Dendrite Nerve Fiber 


Introduction 


In the prior modules on static and dynamic electrical systems, we analyzed basic, hypothetical one-branch nerve 
fibers using a modeling methodology we dubbed the Strang Quartet. You may be asking yourself whether this 
method is stout enough to handle the real fiber of our minds. Indeed, can we use our tools in a real-world setting? 


Presentation 


An Actual Nerve Fiber 


> koe =P a YU SL Oe 


4 oe 


eS Yi 
Le tS se - 


= 


: ~ S - : 
be otreoh Gate en | 
7, Sr pete S aa :- 

Ae TSS a: SAS 


A pyramidal neuron from the CA3 region of a 
rat's hippocampus, scanned at (FIX ME) X 
magnification. 


To answer your question, the above is a rendering of a neuron from a rat's hippocampus. The tools we have refined 
will enable us to model the electrical properties of a dendrite leaving the neuron's cell body. A three-branch model 
of such a dendrite, traced out with painstaking accuracy, appears in the diagram below. 

3-branch Dendrite Model 


Ey Ee Ey 
vd fon Rn Sn { Rin Sn 
¥ 12] { Yas] { ¥. { 
Rj Yial ¥, wv ae 
¥3 ¥e Ys Yio 
Ya ¥is Ry Yue Ry 
R; R; ¥22 R, 22 Ry 
Ra’ Ri a 
Son Ry Ry i 
» 05) } | Sn » So 5 Y¥26 
¥ E, Va ¥| | Ye Sa vl Sn | So 
En En ¥al Yad Ye 
En Ey n 


Multi-compartment electrical model of a rendered dendrite fiber. 


Our multi-compartment model reveals a 3 branch, 10 node, 27 edge structure to the fiber. Note that we have 
included the Nernst potentials, the nervous impulse as a current source, and the additional leftmost edges depicting 
stimulus current shunted by the cell body. 


We will continue using our previous notation, namely: R; and R,,, denoting cell body and membrane resistances, 
respectively; x representing the vector of potentials x. . .z19, and y denoting the vector of currents yj. . .y27. 


Using the typical value for a cell's membrane capacitance, 
c= 1(uF/cm’), 


we derive (see variable conventions): 


Capacitance of a Single Compartment 
l 
C= ana-e 


This capacitance is modeled in parallel with the cell's membrane resistance. Additionally, letting A, denote the 
cell body's surface area, we recall that its capacitance and resistance are 


Capacitance of cell body 
Co a Ape 


Resistance of cell body 
Ra = Acbpm- 


Applying the Strang Quartet 


Step (S1')--Voltage Drops 


Let's begin filling out the Strang Quartet. For Step (S1'), we first observe the voltage drops in the figure. Since 
there are a whopping 27 of them, we include only the first six, which are slightly more than we need to cover all 
variations in the set: 


ej = 71 
€9 = 24, — En 
€3 = f1 — £2 

e4=— 7% 
€5 = 22 — En 
€g5 = 19-23... 
€27 = t10 — En 


In matrix for, letting b denote the vector of batteries, 


onroo7nrcceonrocnrcm;UOUlUvtWrCcOoOmUCcCUOCUmvUMCOCUcOUMWrTWCOWCUOUcUrTWCO CO CUcLUTC OC CO Clr 


(—Em) 


b= 


e—b—Az where 


and 


-l1 0 0 0 0 0 0 0 0 0 
-l1 0 0 0 0 0 0 0 0 0 
-1 1 0 0 0 0 0 0 0 0 
0 -1 0 0 0 0 0 0 0 0 
0 -1 0 0 0 0 0 0 0 0 
0 -1 1 0 0 0 0 0 0 0 
0 0 -1 0 0 0 0 0 0 0 
0 0 -1 0 0 0 0 0 0 0 
0 O -!1 1 0 0 0 0 0 0 
0 0 O 1 -1 0 0 0 0 0 
0 0 0 0 -!1 0 0 0 0 0 
0 0 0 0 -!1 0 0 0 0 0 
0 0 O0 0 -!1 1 0 0 0 0 
A= 0 0 0 0 0 -1 0 0 0 #0 
0 0 0 0 0 -1 0 0 0 0 
0 0 0 0 0 -1 1 0 0 0 
0 0 0 0 0 0 -1 0 0 0 
0 0 0 0 0 0 -1 0 0 0 
0 0 O -1 0 0 0 1 +0 0 
0 0 0 0 0 0 0 -1 0 0 
0 0 0 0 0 0 0 -1 0 0 
0» “Oe “Oo 0° 00 0 0.. = 1. 1. 0 
0 0 0 0 0 0 0 0 -1 0 
0 0 0 0 0 0 0 0 -1 0 
0 0 0 0 0 0 0 0 -1 #1 
0 0 0 0 0 0 0 0 0 -!1 
0 0 0 0 0 0 0 0 0 -!1 


Although our adjacency matrix A is appreciably larger than our previous examples, we have captured the same 
phenomena as before. 


Applying (S2): Ohm's Law Augmented with Voltage-Current Law for Capacitors 


Now, recalling Ohm's Law and remembering that the current through a capacitor varies proportionately with the 
time rate of change of the potential across it, we assemble our vector of currents. As before, we list only enough of 
the 27 currents to fully characterize the set: 


In matrix terms, this compiles to 


y=Ge+C 


where 


Equation: 


Conductance matrix 


oo ° 
oo os 
—) =| =) 
oo ° 


and 


Equation: 


Capacitance matrix 


Cy 90 000 0000 000000000 0 00 0 00 0 


000 000 0000 000 0 00 000 0 00 0 00 0 


000 000 0000 000 0 00 000 0 00 0 00 0 


Cn 00 0000 000000000 0 00 000 0 
000 000 0000 000 0 00 000 0 00 0 00 0 


0 0 0 


000 000 0000 000 0 00 000 0 00 0 00 0 


000 000 000 0 00 000 0 00 0 


Cm 


000 000 0000 000 0 00 000 0 00 0 00 0 


0 00 0 00 


000 000 0000 000 0 00 000 0 00 0 00 0 


000 000 0000 000 0 00 000 0 00 0 00 0 


Cn 00 000 000 000 0 00 0 


000 000 0000 000 0 00 000 0 00 0 00 0 


000 000 0 00 0 


000 000 0000 000 0 00 0 00 0 00 0 00 0 


Cn 00 000 000 0 00 0 


000 000 0000 000 0 00 000 0 00 0 00 0 


000 000 0000 0 00 


C= 


000 000 0000 000 0 00 0 00 0 00 0 00 0 


00 000 0 00 0 


Cm 


000 000 0000 000 0 00 000 0 00 0 00 0 


000 000 0000 000 0 0 0 


000 000 0000 000 0 00 000 0 00 0 00 0 


00 0 00 0 


Cm 


000 000 0000 000 0 00 000 0 00 0 00 0 


000 000 0000 000 0 00 0 00 


000 000 0000 000 0 00 000 0 00 0 00 0 


0 


Cm 0 0 


000 000 0000 000 0 00 000 000 0 00 0 


000 000 0000 000 0 00 0 00 0 00 


000 000 0000 000 0 00 000 0 00 0 00 0 


Cm 


000 000 0000 000 0 00 0 00 0 00 0 00 
000 000 0000 000 0 00 000 0 00 0 00 0 


Step (S3): Applying Kirchoff's Law 


Our next step is to write out the equations for Kirchoff's Current Law. We see: 


Io — Y1 — Yo — y3 = 0 


Y3 — Ys — Ys — Ye = O 


Yo — ¥7 — ys — yo = O 


0 


yo — Y10 — Yig 


0 
0 


Yio — Y11 — Yi2 — Y13 


Yio — Y11 — Yi2 — Y13 


Yi3 — Ya — Ys — Ye = O 
yie — Yi7 — Yis = O 
Yi9 — Y20 — Y21 — Yr2 = O 
Yo2 — Y23 — Y2a — Yrs = O 
Yy25 — Y2 — yor = 0 
Since the B coefficient matrix we'd form here is equal to A’, we can say in matrix terms: 
A’y=—f 


where the vector f is composed of f; = zg and fo..27 = 0. 


Step (S4): Stirring the Ingredients Together 


Step (S4) directs us to assemble our previous toils together into a final equation, which we will then endeavor to 
solve. Using the process derived in the dynamic Strang module, we arrive at the equation 
Equation: 


Arca + ATGAw = ATGD+ f+ aro, 


which is the general form for RC circuit potential equations. As we have mentioned, this equation presumes 
knowledge of the initial value of each of the potentials, 2(0) = X. 


Observing our circuit, and letting a = Goo, we calculate the necessary quantities to fill out [link]'s pieces (for 


these calculations, see dendrite.m): 


ATCA= 


oo o0oco.0UW0mlmlUcCOUcOlhUhO 

SSeS ooo ge Ps 
coo e ee ee 6 
oo oo 0c.o0cdcUmUcOlhlUcOlhUhO 
Secee]ee.c Se 
S680 886 60.6 
e206 22 6-6 ee 
ePpecoocoos 
Sep eooooe cs 


CeCe. -G; 0 0 0 0 0 0 0 
-—G; 2G;4+Gn  —-G; 0 0 0 0 0 0 
0 -G; 2G;+Gn —-G; 0 0 0 0 0 
0 0 =@ 3G; =G; 0 0 =G 0 
7 0 0 0 =¢, 24: +c,, =G 0 0 0 
ios 0 0 0 0 =G, %64C, =C; 0 0 
0 0 0 0 0 =; Gri, 0 0 
0 0 0 =G; 0 0 0 2C.4C, =G 
0 0 0 0 0 0 0 =G, 2G)+4 
0 0 0 0 0 0 0 0 —G 

Go 

Gm 

Gm 

0 

T Gm 

A’Gb = Em G., 

Gn 

Gm 

Gm 

Gin 

7 db 
A Cae = 0, 


and an initial (rest) potential of 


Applying the Backward-Euler Method 


Since our system is so large, the Backward-Euler method is the best path to a solution. Looking at the matrix 
ACA, we observe that it is singular and therefore non-invertible. This singularity arises from the node connecting 
the three branches of the fiber and prevents us from using the simple equation x’ = Bz + g, we used in earlier 
Backward-Euler-ings. However, we will see that a modest generalization to our previous form yields [link]: 
Equation: 


D2a' = Ex+g 


capturing the form of our system and allowing us to solve for a(t). We manipulate [link] as follows: 
Dz' = Ex +g 


#(t) — Z(t — dt) 
dt 


(D — Edt)z(t) = Dz(t — dt) + gdt 


D 


= Ex(t)+g 


%(t) = (D — Edt)‘ (Dz(t — dt) + gdt), 
where in our case 
D=ACA, 
E = — (A’GA), and 
g=A'Gb+ f. 
This method is implemented in dendrite.m with typical cell dimensions and resistivity properties, yielding the 


following graph of potentials. 
Graph of Dendrite Potentials 
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Complex Numbers, Vectors and Matrices 


Complex Numbers 


A complex number is simply a pair of real numbers. In order to stress however that the two 
arithmetics differ we separate the two real pieces by the symbol 7. More precisely, each complex 
number, z, may be uniquely expressed by the combination x + zy, where z and y are real and 2 


denotes / —1. We call zx the real part and y the imaginary part of z. We now summarize the 


main rules of complex arithmetic. 
If 21 = 2, + ty, and zo = x2 + tye then 


Complex Addition 
4y+2=%44+ 224 i(y1 yo) 


Complex Multiplication 


ZZ = (x1 + tyy) (Ho + ty) = 21x_ — Yryo + 1 (H1yo + L2y1) 


Complex Conjugation 
z= 21 —y1 


Complex Division 
cn r1tot+yiyoti(xoyi—x1y2) 


22 oa x74 Yo” 


Magnitude of a Complex Number 


|z1| = Vu = Vt? + yi? 


Polar Representation 


In addition to the Cartesian representation z = x + zy one also has the polar form 


Equation: 


z = |z| (cos(@) +z sin(6)) 


where 0 = arctan(+). 


This form is especially convenient with regards to multiplication. More precisely, 


Equation: 


z1Z2 = |z1||Z2| (cos(@,) cos(@2) — sin(9;) sin(@2) + 7 (cos(0,) sin(@:) + sin(6,) cos(62))) 
= |z1||zZs| (cos(@, + 62) + ¢sin(6, + 62)) 


As a result: 
Equation: 


z” = (|z|)" (cos(n6) + 7 sin(n@)) 


Complex Vectors and Matrices 


A complex vector (matrix) is simply a vector (matrix) of complex numbers. Vector and matrix 
addition proceed, as in the real case, from elementwise addition. The dot or inner product of two 
complex vectors requires, however, a little modification. This is evident when we try to use the 
old notion to define the length of a complex vector. To wit, note that if: 


then 


Now length should measure the distance from a point to the origin and should only be zero for 
the zero vector. The fix, as you have probably guessed, is to sum the squares of the magnitudes 
of the components of z. This is accomplished by simply conjugating one of the vectors. Namely, 
we define the length of a complex vector via: 

Equation: 


(2) = V Fz 


In the example above this produces 


(+a)? + (1a)? = V4 =2 
As each real number is the conjugate of itself, this new definition subsumes its real counterpart. 


The notion of magnitude also gives us a way to define limits and hence will permit us to 


introduce complex calculus. We say that the sequence of complex numbers, 
1 


Zn N= 2 , converges to the complex number Zo and write 


Zn —> 2 
Or 


zy = limit. z, 
n—-> oo 


when, presented with any € > 0 one can produce an integer N for which |z,, — zo| < € when 
7 \n 
n > N. As an example, we note that (4) — 0. 


Examples 


Example: 

As an example both of a complex matrix and some of the rules of complex arithmetic, let us 
examine the following matrix: 

Equation: 


y 
II 
oa 


Let us attempt to find FF’. One option is simply to multiply the two matrices by brute force, 
but this particular matrix has some remarkable qualities that make the job significantly easier. 
Specifically, we can note that every element not on the diagonal of the resultant matrix is equal 
to 0. Furthermore, each element on the diagonal is 4. Hence, we quickly arrive at the matrix 
Equation: 


EE = 


S o& @& fs 
oOo Fk OC 
aS iS S&S © 
Lr Oo CO O&O 


= w% 


This final observation, that this matrix multiplied by its transpose yields a constant times the 
identity matrix, is indeed remarkable. This particular matrix is an example of a Fourier matrix, 
and enjoys a number of interesting properties. The property outlined above can be generalized 
for any F;,, where F' refers to a Fourier matrix with n rows and columns: 

Equation: 


FF, =nl 


Complex Functions 


Complex Functions 


A complex function is merely a rule for assigning certain complex numbers 
to other complex numbers. The simplest (nonconstant) assignment is the 
identity function f(z) = z. Perhaps the next simplest function assigns to 
each number its square, i.e., f(z) = z*. As we decomposed the argument 
of f, namely z, into its real and imaginary parts, we shall also find it 
convenient to partition the value of f, z? in this case, into its real and 
imaginary parts. In general, we write 

Equation: 


f(z + ty) = u(z, y) + iv(z, y) 


where wu and v are both real-valued functions of two real variables. In the 
case that f(z) = z we find 


u(z,y) = 2? -y" 


and 
v(x, y) = 2ary 


With the tools of complex numbers, we may produce complex polynomials 
Equation: 


1 


f(z) = 2” + em-12" +... +e1z+ 9 


We say that such an f is order m. We shall often find it convenient to 
represent polynomials as the product of their factors, namely 
Equation: 


f(z) = (2-1) ™ (z= Ag)®.. «(z= Ap)” 


Each A; is a root of degree d;. Here h is the number of distinct roots of f. 
We call A; a simple root when d; = 1. We can observe the appearance of 
ratios of polynomials or so called rational functions. Suppose 


in rational, that f is of order at most m — 1 while g is of order m with the 
simple roots {A1,..-, Am}. It should come as no surprise that such g should 
admit a Partial Fraction Expansion 

Equation: 


One uncovers the q; by first multiplying each side by z — A; and then 
letting z tend to A;. For example, if 
Equation: 


1 _ ny q_ 
z2+1 z+i z—i 


then multiplying each side by z + 2 produces 


Equation: 
go (z+1 
Z=4 Z=% 
Now, in order to isolate qi it is clear that we should set z = —12. So doing 
we find that gq; = 3. In order to find gz we multiply [link] by z — 7. and 
then set z = 2. So doing we find gz = 3, and so 


Equation: 


a 2 
1 2 2 


2+4 gti et 
. Returning to the general case, we encode the above in the simple formula 
Equation: 

qj = limit (z—Aj;)q(z) 


22M; 


You should be able to use this to confirm that 
Equation: 


2 1/2 1/2 
= ra ; 
27 z+i z-i 


Recall that the transfer function we met in The Laplace Transform module 
was in fact a matrix of rational functions. Now, the partial fraction 
expansion of a matrix of rational functions is simply the matrix of partial 
fraction expansions of each of its elements. This is easier done than said. 
For example, the transfer function of 


is 
Equation: 


(zl — B) = Pal (, ) 


The first line comes form either Gauss-Jordan by hand or via the symbolic 
toolbox in Matlab. More importantly, the second line is simply an 
amalgamation of [link] and [link]. Complex matrices have finally entered 
the picture. We shall devote all of Chapter 10 to uncovering the remarkable 
properties enjoyed by the matrices that appear in the partial fraction 
expansion of (zJ — B) ~' Have you noticed that, in our example, the two 
matrices are each projections, and they sum to J, and that their product is 0? 
Could this be an accident? 


In The Laplace Transform module we were confronted with the complex 
exponential. By analogy to the real exponential we define 
Equation: 


and find that 
Equation: 


e° 


| 
— 
+ 
SN 
+ 
+ 
+ 
at 
+ 


cos(9) + 7 sin(6) 


With this observation, the polar form is now simply z = |z\e’”. 


One may just as easily verify that 


i0 1 ,(-i)0 
cos(9) = ail le 
2 
and 
id —,(--2)0 
sin(9) = “——< 


These suggest the definitions, for complex z, of 


Equation: 
ez te e(t-2) 
cos(z) = 
2 
and 
Equation: 
iz _ p(-i)z 
sin(z) = 
21 


As in the real case the exponential enjoys the property that 


erie —_ pre 2 


and in particular 
Equation: 


erty = etely 


e” cos(y) + te” sin(y) 


Finally, the inverse of the complex exponential is the complex logarithm, 
Equation: 


for z = |z\e*’. One finds that In(—1 + 7%) = in( v2) + 43r, 


Complex Differentiation 


Complex Differentiation 


The complex f is said to be differentiable at zo if 


exists, by which we mean that 


converges to the same value for every sequence {z,,} that converges to zo. In this case we naturally 
call the limit - f (20) 


To illustrate the concept of 'for every' mentioned above, we utilize the following picture. We assume 
the point zo is differentiable, which means that any conceivable sequence is going to converge to Zo. 
We outline three sequences in the picture: real numbers, imaginary numbers, and a spiral pattern of 
both. 

Sequences Approaching A Point In The Complex Plane 


The green is real, the blue is imaginary, and the 
red is the spiral. 


Examples 


Example: 
The derivative of z? is 2z. 


Equation: 
. . . . i +z9) 
limit 2-200 = Jimit “20)te0) 
Z—>Zz0 ay Z—>20 aa) 
= 220 
Example: 
The exponential is its own derivative. 
Equation: 
Z__@20 A z-29 
limit <=" = e*limit =" 
Zz ~m~0 Zz ~—r0 
= elimit 1°, Ga 
elimi i 
pas Se (n+1)! 
= e7 
Example: 


The real part of z is not a differentiable function of z. 
We show that the limit depends on the angle of approach. First, when z,, —> zo ona line parallel to the 


real axis; €:S.5.2,, — @o + + 1y9, we find 
Equation: 


to + + — 29 


limit i - —=1 
PIES Loe = ot ty = Lo 1 10 


while if z, — Zo in the imaginary direction, e.g., 2, = Zo +2 (yo + ae then 
Equation: 


Conclusion 


This last example suggests that when f is differentiable a simple relationship must bind its partial 
derivatives in x and y. 
Partial Derivative Relationship 


If f is differentiable at zo then <4 f(zo) = Ofte) (i Silay) ) 


With z= x+ typ, 


Equation: 
d — Jimie F=f) 
ae f(zo) = limit : 
= aie f(x+tyo)—f(xo+tyo) 
LL t—X0 
_ OF (Zo) 
_ Ox 
Alternatively, when z = zo + 2y then 
Equation: 
da — diene Fl)=fl20) 
dz f(z) = limit Z—29 
= limit f(xo+iy)—f(xo+tyo) 
y— Yo i(y—Yo) 


Cauchy-Reimann Equations 


In terms of the real and imaginary parts of f this result brings the Cauchy-Riemann equations. 
Equation: 


A eeegcts 
Ox Oy 
and 
Equation: 
OM 2.08 
Ox ~—S Oy 


Regarding the converse proposition we note that when f has continuous partial derivatives in region 
obeying the Cauchy-Reimann equations then f is in fact differentiable in the region. 


We remark that with no more energy than that expended on their real cousins one may uncover the 
rules for differentiating complex sums, products, quotients, and compositions. 


As one important application of the derivative let us attempt to expand in partial fractions a rational 
function whose denominator has a root with degree larger than one. As a warm-up let us try to find q11 
and q1,2 in the expression 


z+2 0, M1 41,2 


(241)? z+1 °° (241)? 


Arguing as above, it seems wise to multiply through by (z + Nig and so arrive at 
Equation: 


z4+2=qyi(z+1)+q2 


On setting z = —1 this gives q1,2 = 1. With qi,2 computed, [link] takes the simple form 
z+1=qii(z+1) andso qi2 = 1 as well. Hence, 


z+2 1 1 


gai 242 Gai? 
This latter step grows more cumbersome for roots of higher degrees. Let us consider 


(z+ 2) 1,1 q1,2 is 91,3 


get el Gai? @tiy 


The first step is still correct: multiply through by the factor at its highest degree, here 3. This leaves us 
with 
Equation: 


(2+ 2)? = qa(zt+1)?+q2(z+1) +03 


Setting z = —1 again produces the last coefficient, here q1,3 = 1. We are left however with one 
equation in two unknowns. Well, not really one equation, for [link] is to hold for all z. We exploit this 
by taking two derivatives, with respect to z, of [link]. This produces 


2(z+2) =2qii(z+1) 4+ M2 


and 2 = qi, The latter of course needs no comment. We derive qi,2 from the former by setting z = —1 
. This example will permit us to derive a simple expression for the partial fraction expansion of the 
general proper rational function g = - where g has h distinct roots {A1,..., An} of respective degrees 
{d,,...,d;,}. We write 

Equation: 


and note, as above, that qj, is the coefficient of (z — d;)“* in the rational function 


rj(z) = a(z)(z — 4)” 


Hence, q;,4 may be computed by setting z = A; in the ratio of the d; — kth derivative of r; to 
(d; — k)!. That is, 


Equation: 


1 d4i-k 
dj.k = limit {(z 


zr, (d; = k)! dzti-* 


As a second example, let us take 


Equation: 
1 0 0 
B= 1 3 0 
1 1 
and compute the @; , matrices in the expansion 
1 
7 0 0 
= 1 1 1 
(I-B)*= TGpyeD Z3 0 = 
1 1 1 


(z-1)°(z—-3) — (z-1)(z-3) —z-1 


The only challenging term is the (3, 1) element. We write 


1 dit 1,2 
+ 


d; 
dj)“ a(2)} 
ie 
z — 
q2,1 


(2-19 (2-3) 2-1 (2-1) 2-3 


It follows from [link] that 


Equation: 
71,1 = ( = 1) 
—1/4 
and 
Equation: 
12 = + 1 
—1/4 
and 
Equation: 


gi = ( sy y"4 
1/4 


It now follows that 
Equation: 


1 0 0 0 00 , 9 0 0 
-1/2 0 0 +——,> 0 00 + 1/2 1 0 
1/4 1721 ©-Y 1200 ” 1/4 1/2 0 


I-B)t= 
io ) gol 


In closing, let us remark that the method of partial fraction expansions has been implemented in 


Matlab. In fact, [link], [link], and [link] all follow from the single command: 
[r,p,k]=residue([0 © 0 1],[1 -5 7 -3]). The first input argument is Matlab-speak for 


the polynomial f(z) = 1 while the second argument corresponds to the denominator 


g(z) = (z-1)? (2-3) = 2-52? 4+ 72-3 


Exercises: Complex Numbers, Vectors, and Functions 


Exercises 


Exercise: 


Problem: Express |e*| in terms of x and/or y. 


Solution: 
Pending completion of assignment. 


Exercise: 


Problem: Confirm that e™*) = z and In(e”) = z. 


Solution: 


Pending completion of assignment. 
Exercise: 


Problem: 


Find the real and imaginary parts of cos(z) and sin(z). Express your 
answers in terms of regular and hyperbolic trigonometric functions. 


Solution: 


Pending completion of assignment. 


Exercise: 


Problem: Show that cos?(z) + sin?(z) = 1 


Solution: 


Pending completion of assignment. 


Exercise: 


Problem: With z” = e”™©) for complex z and w compute Vi 


Solution: 


Pending completion of assignment. 
Exercise: 


Problem: 


Verify that cos(z) and sin(z) satisfy the Cauchy-Riemann equations 
and use the proposition to evaluate their derivatives. 


Solution: 


Pending completion of assignment. 
Exercise: 
Problem: 


Submit a Matlab diary documenting your use of residue in the partial 
fraction expansion of the transfer function of 


2 0 @ 
B= -1 4 0 
0 =-l1 2 


Solution: 


Pending completion of assignment. 


Cauchy's Theorem 


Introduction 


Our main goal is a better understanding of the partial fraction expansion of a given transfer function. With respect 
to the example that closed the discussion of complex differentiation, see this equation - In this equation, we found 


1 1 1 
Pi D4 P. 
Pe ia eo? n” gael 


(zI — B)* = 


where the P; and D; enjoy the amazing properties 


1. Equation: 
BP, = PB 
= AP,+D, 
and 
BP: = PoB = 2P» 
2. Equation: 
Pi+P,=I1 
Pi? =P, 
P,? = Py 
and 
De =0 
3. Equation: 
PD, = DiPi 
dD, 
and 


P,D, = D, Py, =0 


In order to show that this always happens, i.e., that it is not a quirk produced by the particular B in this equation 
we require a few additional tools from the theory of complex variables. In particular, we need the fact that partial 
fraction expansions may be carried out through complex integration. 


Integration of Complex Functions Over Complex Curves 


We shall be integrating complex functions over complex curves. Such a curve is parameterized by one complex 
valued or, equivalently, two real valued, function(s) of a real parameter (typically denoted by t). More precisely, 


C = {z(t) = a(t) +iy(t)|a<t <b} 
For example, if x(t) = y(t) = t while a = 0 and b = 1, then C is the line segment joining 0 + 70 to 1 +2. 


We now define 


[rears [reat 


For example, if C = {¢ + it|0 < t < 1} as above and f(z) = z then 


1 L 
feaz=f +i a+adt=f t—t+i2tidt=i 
0 0 


while if C is the unit circle {e”|0 < t < 27} then 


Qn Qn Qn 
feaz-f ctitat =i [ Matai f cos(2t) + isin(2t) dt = 0 
0 0 0 


Remaining with the unit circle but now integrating f(z) = 4+ we find 


Qn 
/ ztdz= | e “et dt = 2mi 
0 


We generalize this calculation to arbitrary (integer) powers over arbitrary circles. More precisely, for integer m 
and fixed complex a we integrate (z — a)" over 


C(a,r) = {a+ re"|0<t < 2n} 


the circle of radius r centered at a. We find 
Equation: 


f(z-a)"dz = fP" (at+re*—a)”rie# dt 


irmtl {* eil(m+lt qt 


2ni if m=—-1 


/ (z—a)"dz=ir™ [F costem + 1)t) +isin((m+1)t)dt= { 


0 otherwise 


When integrating more general functions it is often convenient to express the integral in terms of its real and 
imaginary parts. More precisely 


[rede [ uey+inley dati f u(e,y) +iv(e,v) dy 


[r@az=[ueryae- foenav+if veyaerif weyay 


/ f(z) dz= / u(x(t), y(t))2'(t) — o(a(t), y(t))y'(t) dt + if u(x(t), y(t))y'(t) + v(x(t), y(t))2"(t) dt 


a 


The second line should invoke memories of: 
Green's Theorem 


If C is a closed curve and M and WN are continuously differentiable real-valued functions on Cin, the region 


enclosed by C, then 
N M 
[act fway=f [5 e dzady 
Ox Oy 


Applying this to the situation above, we find, so long as C is closed, that 


| faz- -[{F 24 acdytif fF + Gr dedy 


At first glance it appears that Green's Theorem only serves to muddy the waters. Recalling the Cauchy-Riemann 
equations however we find that each of these double integrals is in fact identically zero! In brief, we have proven: 
Cauchy's Theorem 


If f is differentiable on and in the closed curve C then f f(z) d z=0. 


Strictly speaking, in order to invoke Green's Theorem we require not only that f be differentiable but that its 
derivative in fact be continuous. This however is simply a limitation of our simple mode of proof; Cauchy's 
Theorem is true as stated. 


This theorem, together with [link], permits us to integrate every proper rational function. More precisely, if 

q= £ where f is a polynomial of degree at most m — 1 and g is an mth degree polynomial with h distinct zeros 
at {A;| 7 = {1,...,h}} with respective multiplicities of {m,| 7 = {1,...,h}} we found that 

Equation: 


h mM; 


q(z) = Dee e ae )F 


g=] k=l 


Observe now that if we choose r; so small that A; is the only zero of g encircled by C; = C(Aj, rj) then by 


Cauchy's Theorem 
m3 1 
g(z)dz= S_ aie ——; dz 
kal (z— Aj) 


In [link] we found that each, save the first, of the integrals under the sum is in fact zero. Hence, 
Equation: 


i q(z) d z= 27ig;1 


With q; in hand, say from this equation or residue, one may view [link] as a means for computing the 
indicated integral. The opposite reading, i.e., that the integral is a convenient means of expressing q;,1, will prove 
just as useful. With that in mind, we note that the remaining residues may be computed as integrals of the product 
of q and the appropriate factor. More precisely, 

Equation: 


f aee- rsh de = aig, 


One may be led to believe that the precision of this result is due to the very special choice of curve and function. 
We shall see ... 


Cauchy's Integral Formula 


The Residue Theorem 


After Cauchy's Theorem egn6 and Cauchy's Theorem egn7 perhaps the most useful consequence 
of Cauchy's Theorem is the 
The Curve Replacement Lemma 


Suppose that C2 is a closed curve that lies inside the region encircled by the closed curve C. If f 
is differentiable in the annular region outside C2 and inside C, then 


[t@az=fteoaz 


With reference to [link] we introduce two vertical segments and define the closed curves 
C’3 = abcda (where the bc arc is clockwise and the da arc is counter-clockwise) and C4 = adcba 
(where the ad arc is counter-clockwise and the cb arc is clockwise). By merely following the 


arrows we learn that 
[r@az=fieazt f teQazs f teaz 


As Cauchy's Theorem implies that the integrals over C’3 and C’4 each vanish, we have our result. 
Curve Replacement Figure 


The Curve Replacement Lemma 


This Lemma says that in order to integrate a function it suffices to integrate it over regions where 
it is singular, i.e. nondifferentiable. 


Let us apply this reasoning to the integral 


lame” 


where C’ encircles both A; and Ag as depicted in [link]. We find that 


lene eee) eee 


Developing the integrand in partial fractions we find 


z yi ; 1 AG / 1 Qrid, 
dz= dz+ dz2= 
ceesnicesn ‘ Ai — A2 z—-r\ ’ Ag — A z—do 7 Ai — Ag 


Similarly, 


/ z F| amir 
z= 
(z—A1) (z— 2) r2 —r, 


Putting things back together we find 
Equation: 


z . r 
J Trj)e-ry) 12% = 2m (x24 og iz) 


= 27 


Image not finished 


Concentrating on 
the poles. 


We may view [link] as a special instance of integrating a rational function around a curve that 
encircles all of the zeros of its denominator. In particular, recalling that cauchy's theorem eqn5 and 
Cauchy's theorem egn6, we find 

Equation: 


| 


fa(z)dz ~ we Cae dz 


. h 


f(z) 


zZ—-a 
of which f is differentiable and a resides. Our Curve Replacement Lemma now permits us to 
claim that 


over some closed curve C’ inside 


To take a slightly more complicated example let us integrate 


ars [az 


z—a 


It appears that one can go no further without specifying f. The alert reader however recognizes 
that in the integral over C(a, 1) is independent of r and so proceeds to let r — 0, in which case 
z—> aand f(z) — f(a). Computing the integral of — along the way we are led to the hope that 


IO) a, = $(a)2ni 
Z= a 
In support of this conclusion we note that 
JO q,- f 0, £0 10) a, = yo f oars f OM a, 
Ze Z=O0 °°. =. 24 yas) 2 


Now the first term is f(a)277 regardless of r while, as r — 0, the integrand of the second term 

approaches 4 (a) and the region of integration approaches the point a. Regarding this second 
term, as the integrand remains bounded as the perimeter of C(a,7) approaches zero the value of 
the integral must itself be zero. This result if typically known as 

Formula 

Cauchy's Integral Formula 


If f is differentiable on and in the closed curve C’ then 
Equation: 


ee ee 


2771 z—a 


dz 


for each a lying inside C. 


The consequences of such a formula run far and deep. We shall delve into only one or two. First, 
we note that, as a does not lie on C, the right hand side is a perfectly smooth function of a. Hence, 
differentiating each side, we find 

Equation: 


for each a lying inside C’. Applying this reasoning n times we arrive at a formula for the n-th 
derivative of f at a, 
Equation: 


d” n! f(z) 
f(a) = dz 
rIO=s/[ 


277% 


for each a lying inside C’. The upshot is that once f is shown to be differentiable it must in fact be 
infinitely differentiable. As a simple extension let us consider 


1 f(z) 


= fe 


Qmi J (z—1)(z— Az)? 


where f is still assumed differentiable on and in C and that C encircles both A; and A». By the 
curve replacement lemma this integral is the sum 


&. fz) oe fl) 7 
wat | (z — A1)(z— Aa)? vee iat | (z — A1)(z— Aa)? ° 


f(z) 


z—Ag 


where A, now lies in only C;. As is well behaved in C; we may use [link] to conclude that 


a ee 


Qi - 


(z— A1)(z— Az)? (Ar — A2)? 
f(z) 


Z—-Xy 


Similarly, as is well behaved in C’2 we may use [link] to conclude that 


= / (z— a 2)? aaae eC ( sas 


These calculations can be read as a concrete instance of 
The Residue Theorem 


a=, 


If g is a polynomial with roots {A,| 7 = {1,..., h}} of degree {d,| 7 = {1,...,h}} andCisa 
closed curve encircling each of the A, and f is differentiable on and in C’ then 


/ 


f(z) z= 2m : res (A, 


where 


res (A;) = to Ce a 


is called the residue of f at A;. 


One of the most important instances of this theorem is the formula for the Inverse Laplace 
Transform. 


The Inverse Laplace Transform: Complex Integration 


The Inverse Laplace Transform 


If q is a rational function with poles {A,| 7 = {1,...,}} then the inverse Laplace transform of gq is 
Equation: 


¥\(q)(t) = — | a(2e* dz 


Qi 


where C is a curve that encloses each of the poles of g. As a result 
Equation: 


(q(t) = Sores (25) 


j=l 


Let us put this lovely formula to the test. We take our examples from discussion of the Laplace Transform and the 
inverse Laplace Transform. Let us first compute the inverse Laplace Transform of 


1 
2) =; 
(z+1) 
According to [link] it is simply the residue of q(z)e* at z = —1, ie., 
a d et 3 
res (—1) = limit, es te 


This closes the circle on the example begun in the discussion of the Laplace Transform and continued in exercise 
one for chapter 6. 


For our next example we recall 


0.19 s?+1.55s+0.27 


L(21(s)) = ; 
(s + 1/6)* (s3 + 1.6555? + 0.49785 + 0.0039) 


from the Inverse Laplace Transform. Using numde, sym2poly and residue, see fib4.m for details, returns 


0.0029 
262.8394 
—474.1929 
ry =  —1.0857 
—9.0930 
—0.3326 
211.3507 


and 


—1.3565 
—0.2885 
—0.1667 
pi= —0.1667 
—0.1667 
—0.1667 
—0.0100 


You will be asked in the exercises to show that this indeed jibes with the 


= -t —329t Vv 73t 
a i(t) = 211.35e% — 0.0554t? + 4.5464t? + 1.085¢ + 474.19 e+ +e 4 262.842cosh ——— + 262. 
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achieved in the Laplace Transform via ilaplace. 


Exercises: Complex Integration 


1. Let us confirm the representation of this Cauchy's Theorem equation in 
the matrix case. More precisely, if 6(z) = (zI — B)~’ is the transfer 
function associated with B then this Cauchy's Theorem equation states 


that 
hd; 
A] Pp “ke 
B(z)= > ; 
j=l k=l (Z— dj)" 
where 
Equation: 
1 k-1 
®;, = —~ | B(z)(z-A;)” dz 
271 


Compute the @; ;, per [link] for the B in this equation from the 
discussion of Complex Differentiation. Confirm that they agree with 
those appearing in this equation from the Complex Differentiation 
discussion. 

2. Use this inverse Laplace Transform equation to compute the inverse 
Laplace transform of aoe 

3. Use the result of the previous exercise to solve, via the Laplace 
transform, the differential equation 


<(a)(t) +2(#) =e*sin(#), 2(0) =0 
Hint: Take the Laplace transform of each side. 

4. Explain how one gets from r, and pj to 2; (t). 

5. Compute, as in Fib4.m, the residues of Y(x2(s)) and #(x3(s)) 
and confirm that they give rise to the x2(t) and x3(t) you derived in 
the discussion of Chapter 1. 


The Eigenvalue Problem 


Introduction 


Harking back to our previous discussion of The Laplace Transform we labeled the 
complex number \ an eigenvalue of B if AJ — B was not invertible. In order to 
find such 4 one has only to find those s for which (sI — B)" is not defined. To 
take a concrete example we note that if 


Equation: 
1 10 
B= 01 0 
00 2 
then 
Equation: 
; (s —1)(s— 2) s—2 0 
(I - By °=——_ (s —1)(s —2) 0 
(s-1)°(s—2) . a0 


and so Ay = 1 and Az = 2 are the two eigenvalues of B. Now, to say that 
A; — B is not invertible is to say that its columns are linearly dependent, or, 
equivalently, that the null space (A; — B) contains more than just the zero 
vector. We call (A,;J — B) the jth eigenspace and call each of its nonzero 
members a jth eigenvector. The dimension of | (A,;J — B) is referred to as the 
geometric multiplicity of \;. With respect to B above, we compute 

(Ail — B) by solving (J — B)zx = 0, ie., 


0 —-1 0 Hi 0 
0 0 0 to = 0 
0 0 1 X3 0 


Clearly 


Arguing along the same lines we also find 
(AI - B)= c(0 0 1)’ ce 


That B is 3x3 but possesses only 2 linearly eigenvectors leads us to speak of B as 
defective. The cause of its defect is most likely the fact that Az is a double pole of 
(sI — B)' In order to flesh out that remark and uncover the missing eigenvector 
we must take a much closer look at the transfer function 


R(s) = (sI — B)* 


In the mathematical literature this quantity is typically referred to as the 
Resolvent of B. 


Eigenvalue Problem: The Transfer Function 


The Transfer Function 


One means by which to come to grips with R(s) is to treat it as the matrix analog 
of the scalar function 
Equation: 


1 
s—b 


This function is a scaled version of the even simpler function 7S This latter 
function satisfies the identity (just multiply across by 1 — z to check it) 
Equation: 

1 e 


—=1 Cs sg 
—— rere a t Z a ae 


for each positive integer n. Furthermore, if |z| < 1 then z” — 0as n — oo and 
so [link] becomes, in the limit, 


al ae 
rier) 


n=0 


the familiar geometric series. Returning to [link] we write 


om [Re 


ate oe: 
sz gn s” s—b 


: 
8 


and hence, so long as |s| > |b| we find, 
1 23 b\" 
s—b 5s = ee 


This same line of reasoning may be applied in the matrix case. That is, 
Equation: 


: B\'/1 B Br-1 Bn : 
(ot - By? =s(1- =) (c+ote + —(sI — B) ) 
S S s” s” 


and hence, so long as |s| >|| B || where || B || is the magnitude of the largest 
element of B, we find 
Equation: 


Although [link] is indeed a formula for the transfer function you may, regarding 
computation, not find it any more attractive than the Gauss-Jordan method. We 
view [link] however as an analytical rather than computational tool. More 
precisely, it facilitates the computation of integrals of R(s). For example, if C,, is 
the circle of radius p centered at the origin and p >|] B || then 

Equation: 


le GI=B) “de = S048" So, ds 
= 2m 


This result is essential to our study of the eigenvalue problem. As are the two 
resolvent identities. Regarding the first we deduce from the simple observation 


(f= BB) (iB) = (f= BY f= B sl +B) By 


that 
Equation: 


R(s2) — R(s1) = (81 — 82) R(s2)R(s1) 


The second identity is simply a rewriting of 
(sI — B)(sI — B)"' =(sI— B)* (sI— B) =I 


namely, 
Equation: 


BR(s) 


R(s)B 
sR(s)—I 


Partial Fraction Expansion of the Transfer Function 


Partial Fraction Expansion of the Transfer Function 


The Gauss-Jordan method informs us that R will be a matrix of rational functions with a common denominator. In 
keeping with the notation of the previous chapters, we assume the denominator to have the h distinct roots, 
{A;|j = {1,..., h}} with associated multiplicities {m,| j = {1,..., h}}. 


Now, assembling the partial fraction expansions of each element of R we arrive at 
Equation: 


where, recalling the equation from Cauchy's Theorem, the matrix R;;, equals the following: 
Equation: 


1 = 
Rik = a | Re@le-Ay ‘dz 


Example: 
Concrete Example 
As we look at this example with respect to the eigenvalue problem eqn1 and eqn2, we find 


100 Oo t @ 0 0 0 
Ri = 010 Ri» = 00 0 and Roi = 00 0 
0 0 0 00 0 0 @ i 


One notes immediately that these matrices enjoy some amazing properties. For example 
Equation: 


Riv? =Rii, Roai?= Rei, RiiRei=0, and Roi7=0 


Below we will now show that this is no accident. As a consequence of [link] and the first resolvent identity, we 
shall find that these results are true in general. 


Rj? = Rj, as seen above in [link]. 


Recall that the C’; appearing in [link] is any circle about A, that neither touches nor encircles any other root. 
Suppose that C; and C;,’ are two such circles and C;’ encloses C;. Now 


1 1 
Ry = 5 f Re)dz= 5 | Raz 


and so 


ig’ (2 aol [Re w)dwdz 
Th) 
R(w 
Ry" = sa a = Pere 
(27) wz 
1 1 
— [ Rw f Se azaw) 
(277) wz 
~ Ont =; | RO) das Red 


We used the first resolvent identity, This Transfer Function eqn, in moving from the second to the third line. In 
moving from the fourth to the fifth we used only 
Equation: 


and 


dz=0 


/ 1 
wz 


The latter integrates to zero because C’; does not encircle w. 


From the definition of orthogonal projections, which states that matrices that equal their squares are projections, 
we adopt the abbreviation 


& = Ry 


With respect to the product P;.P;, for 7 # k, the calculation runs along the same lines. The difference comes in 
[link] where, as C; lies completely outside of C’;, both integrals are zero. Hence, 


If 7 A k then P,P, = 0. 
Along the same lines we define 
D; = Rj2 
and prove 
If1<k<m,;—1 then D;* = Rjyy. Dj =0. 


For k and / greater than or equal to one, 


a z\(z—r;)*dz w)(w —A;)' dw 
RyneaRuea = oe | RO ride f R(w)w—A,)'d 


Rp riRji1 = 


w)(z— dj) *(w - dy)! dwdz 


ee R(2) ~ Rw) iy awa e 
Ripe Rj = (nip? I = (z—Aj;)"(w—A;) dwd 


RieiRyns = f R@(e-A)" [ GW aware [ reoyo—ayt f 22" 


1 
Rye Rjis = Oni / R(z)(z = oe bal dz= Rye 


because 
Equation: 
w—,) 
/ iia dw = 2ni(z — A,)! 
wz 
and 


[2 emo 


Wz 


With k = 1 = 1 we have shown Rj2” = Rj3, i.e., Dj” = Rj3. Similarly, with k = 1 and = 2 we find 
Rj2Rj3 = Rj, ie, Dj? = Rj. Continuing in this fashion we find Rj ~Rj4041 = Rj pro = J, or 
ae = Rjx+2. Finally, at k = m; — 1 this becomes 


1 
Dj” = Rymyi = 5 | Rl) —Aj)™ dz=0 


27% 
by Cauchy's Theorem. 


With this we now have the sought after expansion 
Equation: 


along with the verification of a number of the properties laid out in Complex Integration eqns 1-3. 


The Spectral Representation 


With just a little bit more work we shall arrive at a similar expansion for B 
itself. We begin by applying the second resolvent identity to P;. More 
precisely, we note that the second resolvent identity implies that 
Equation: 


== Jo, #R(z) —Idz 


P23 / zR(z) dz 
é 


PBs | R(z) (z— 2;) dz, [ ee 
Cj Cj 
PPa Dae 


Summing this over 7 we find 
Equation: 


We can go one step further, namely the evaluation of the first sum. This 
stems from the eqn in the discussion of the transfer function where we 
integrated R(s) over a circle C’, where p >|| B ||. The connection to the P; 
is made by the residue theorem. More precisely, 


i 


p 


h 
R(z)dz=2ni)- P, 
j=l 


Comparing this to the eqn from the discussion of the transfer function we 
find 


Equation: 


j=l 
and so [link] takes the form 
Equation: 
h h 
B=S°A;P;/+S_D; 
j=1 j=l 


It is this formula that we refer to as the Spectral Representation of B. To 
the numerous connections between the P; and D, we wish to add one more. 
We first write [link] as 


(B—A,I)P; =D; 


and then raise each side to the m,; power. As Ps = P; and D = 0 we 
find 
Equation: 


(B—,I)™P; =0 


For this reason we call the range of P; the 7th generalized eigenspace, call 
each of its nonzero members a 7th generalized eigenvector and refer to the 
dimension of (P,) as the algebraic multiplicity of A;. With regard to the 
first example from the discussion of the eigenvalue problem, we note that 
although it has only two linearly independent eigenvectors the span of the 
associated generalized eigenspaces indeed fills out *. One may view this 
as a consequence of P; + P) = I, or, perhaps more concretely, as 
appending the generalized first eigenvector (0 1 0)" to the original two 


eigenvectors (1 0 0)’ and(0 0 1)’. Instill other words, the algebraic 


multiplicities sum to the ambient dimension (here 3), while the sum of 
geometric multiplicities falls short (here 2). 


The Eigenvalue Problem: Examples 


We take a look back at our previous examples in light of the results of two 
previous sections The Spectral Representation and The Partial Fraction 
Expansion of the Transfer Function. With respect to the rotation matrix 


e- (4 o) 


we recall, see Cauchy's Theorem egné6, that 


RO) = a5 « , 


Equation: 


and so 


17? = 12. = 
p= apr ap=(% Ae 4 
5 1/2 = 12 


From m,; = mz = 1 it follows that @(P,) and &(P2) are actual (as 
opposed to generalized) eigenspaces. These column spaces are easily 
determined. In particular, @(P;) is the span of 


~=() 


while &(P2) is the span of 


i C) 


To recapitulate, from partial fraction expansion one can read off the 
projections from which one can read off the eigenvectors. The reverse 
direction, producing projections from eigenvectors, is equally worthwhile. 
We laid the groundwork for this step in the discussion of Least Squares. In 


P,= e1(e1"e1) ee and P, =e» (e2"e2) ee 
As e; e, = e;" €; = 0 these formulas can not possibly be correct. 
Returning to the Least Squares discussion we realize that it was, perhaps 
implicitly, assumed that all quantities were real. At root is the notion of the 
length of a complex vector. It is not the square root of the sum of squares of 
its components but rather the square root of the sum of squares of the 
magnitudes of its components. That is, recalling that the magnitude of a 


complex quantity z is V 22, 


({| e1 ||)” A ex"e1 rather (| e1 ||)" 4 Ger 


Yes, we have had this discussion before, recall complex numbers, vectors, 
and matrices. The upshot of all of this is that, when dealing with complex 
vectors and matrices, one should conjugate before every transpose. Matlab 
(of course) does this automatically, i.e., the ' symbol conjugates and 
transposes simultaneously. We use x" to denote ‘conjugate transpose’, i.e., 
at = gt 
All this suggests that the desired projections are more likely 
Equation: 


P, = Ee] (evei) e™ and P, = €2 (e2"e2) és 


Please check that [link] indeed jives with [link]. 


The Eigenvalue Problem: Exercises 


Exercises 


1. Argue as in Proposition 1 in the discussion of the partial fraction 
expansion of the transfer function that if 7 4 k then 
D,P, = P;D, = 0. 

2. Argue from this equation from the discussion of the Spectral 
Representation that D;P; = P;D; = Dj. 

3. The two previous exercises come in very handy when computing 
powers of matrices. For example, suppose B is 4-by-4, that h = 2 and 
my, = Mz = 2. Use the spectral representation of B together with the 
first two exercises to arrive at simple formulas for B? and B?. 

4. Compute the spectral representation of the circulant matrix 


Carefully label all eigenvalues, eigenprojections and eigenvectors. 


The Spectral Representation of a Symmetric Matrix 


Introduction 
Our goal is to show that if B is symmetric then 


e each A; is real, 
e each P; is symmetric and 
e each D; vanishes. 


Let us begin with an example. 


Example: 
The transfer function of 


is 


WG ile ie , 
Re ee se 
We Sue aye Ye 


1/33 
yay 
1 31/3 


and so indeed each of the bullets holds true. With each of the D; falling by 
the wayside you may also expect that the respective geometric and 
algebraic multiplicities coincide. 


The Spectral Representation 


We have amassed anecdotal evidence in support of the claim that each D; 
in the spectral representation 
Equation: 


is the zero matrix when B is symmetric, i.e., when B = BT, or, more 
generally, when B = B" where BY = B® Matrices for which B = B# 
are called Hermitian. Of course real symmetric matrices are Hermitian. 


Taking the conjugate transpose throughout [link] we find 
Equation: 


That is, the A; are the eigenvalues of B#* with corresponding projections 
P;# and nilpotents D;# Hence, if B = B™, we find on equating terms that 


Az = Aj 

PP 
and 

D.=pD# 


The former states that the eigenvalues of an Hermitian matrix are real. Our 
main concern however is with the consequences of the latter. To wit, notice 
that for arbitrary x, 


| D.mi1¢ | 2 _ Pa D.mi-1 Hy mj-1,, 
J Fi) a 
| D;"i1¢ | 2 - zip .mi—1 pF) .mj—1,, 
J J J 
DD. 1¢ | 2 — zi D .Mj—2 J) Mj -p 
(|| D; | j j 
(|| D;""*2 ||)" =0 


As D Pima = 0 for every x it follows (recall this previous exercise) that 
Dy “l_Qg, Continuing in this fashion we find D;"” ~2 — 0 and SO, 
eventually, D; = 0. If, in addition, B is real then as the eigenvalues are real 
and all the D; vanish, the P; must also be real. We have now established 


If B is real and symmetric then 
Equation: 


where the A; are real and the P; are real orthogonal projections that sum to 
the identity and whose pairwise products vanish. 


One indication that things are simpler when using the spectral 
representation is 
Equation: 


h 
pio = S- dj PP; 
j=l 


As this holds for all powers it even holds for power series. As a result, 


h 
e? = ) iP, 
j=l 


It is also extremely useful in attempting to solve 
Be=b 


for x. Replacing B by its spectral representation and b by Jb or, more to the 
point by >); P;b we find 


WE 
> 
ae 

8 

| 
M 
by 

ao 


j= Via 


Multiplying through by P; gives A; P,z2 = P,b or Pix = oe Multiplying 
through by the subsequent P;,'s gives P;x = Hence, 
Equation: 


h 
t= dij Pie 
h 
=v 


We clearly run in to trouble when one of the eigenvalues vanishes. This, of 
course, is to be expected. For a zero eigenvalue indicates a nontrivial null 
space which signifies dependencies in the columns of B and hence the lack 
of a unique solution to Bx = b. 


Another way in which [link] may be viewed is to note that, when B is 
symmetric, this previous equation takes the form 


= 
a 


(II-B) *=S- 


Now if 0 is not an eigenvalue we may set z = 0 in the above and arrive at 
Equation: 


| 
=) 


j=l °I 


Hence, the solution to Bx = bis 


h 
_ 1 
j=l Jj 


as in [link]. With [link] we have finally reached a point where we can begin 
to define an inverse even for matrices with dependent columns, i.e., with a 
zero eigenvalue. We simply exclude the offending term in [link]. Supposing 
that A;, = 0 we define the pseudo-inverse of B to be 


Let us now see whether it is deserving of its name. More precisely, when 
b € &(B) we would expect that c = Bb indeed satisfies Br = b. Well 


> h-1 1 h-1 
BB*b=BY_ 5, Pib = 
j=l “J 


It remains to argue that the latter sum really is b. We know that 


h 
Vb,bER(B): b= > ° Pb 


j=l 


The latter informs us that b | N Toaee As B = B’, we have, in fact, that 
b | N(B). As Py, is nothing but orthogonal projection onto N(B) it 
follows that P;,b = 0 and so B(B*b) = 6, that is, r = B*b is a solution to 


Bz = b. The representation [link] is unarguably terse and in fact is often 
written out in terms of individual eigenvectors. Let us see how this is done. 
Note that if ¢ € R(P;) then z = P,x and so, 


h 
be= BP x = SC A;P)Pix = Ay P 1x = Aix 
j=l 


i.e., z is an eigenvector of B associated with A. Similarly, every (nonzero) 
vector in S4(P;) is an eigenvector of B associated with ,. 


Next let us demonstrate that each element of S%(P;) is orthogonal to each 
element of R(P;,) when j # k. If « € K(P;) and y € R(P;) then 


aly — (P,x)" Pry = a’ P;Pry — 0 


With this we note that if 1e Cea pee 2 pay constitutes a basis for 
$(P,) then in fact the union of such bases, 


{tjp| (1 Sj <h) \ As p<nj)} 
forms a linearly independent set. Notice now that this set has 
h 
ie 
j=l 


elements. That these dimensions indeed sum to the ambient dimension, n, 
follows directly from the fact that the underlying P; sum to the n-by-n 
identity matrix. We have just proven 


If B is real and symmetric and n-by-n, then B has a set of n linearly 
independent eigenvectors. 


Getting back to a more concrete version of [link] we now assemble matrices 
from the individual bases 


i {e515 Di Q ee) 2 syn; } 


and note, once again, that P; = £; (E;" E;) a oF and so 


h 
B=)_ jE; (E;7E;) 'E;? 
j=l 


I understand that you may feel a little overwhelmed with this formula. If we 
work a bit harder we can remove the presence of the annoying inverse. 
What I mean is that it is possible to choose a basis for each #(P;) for 
which each of the corresponding F;; satisfy E;’ EF; = I As this 
construction is fairly general let us devote a separate section to it (see 
Gram-Schmidt Orthogonalization). 


Gram-Schmidt Orthogonalization 

Suppose that MV is an m-dimensional subspace with basis 
{Hisesag 

We transform this into an orthonormal basis 
{q lye esy Qm} 

for M via 


1, Set yy = a7 and:q; = Th 


2. Y2 = £2 minus the projection of x2 onto the line spanned by qj. That 
is 


-1 
Y2 = %2— 1 (q1° 41) qi 2 — a gig x2 


Set q2 = La and Q2 = [q1, 42]. 
3. y3 = £3 minus the projection of x3 onto the plane spanned by q; and 
q2. That is 


= 
y3 = 23 — Qe (Q2*Q2) Qo’ a3 = 3 = qq x3 


Set q3 = Tal and Q3 = {q1, 92, 93}. Continue in this fashion through 


step (m). 


¢ (Mm) Ym = Lm Minus its projection onto the subspace spanned by the 
columns of Q,,_1. That is 


m—1 
Ym = Lm — Omak (Qm-1° Qm-1) Oni tte, = 934; 2m 
ql 


Set dm = i To take a simple example, let us orthogonalize the 
3 


following basis for °, 


lqa=y1= 7}. " 
2. yo = 22—-q191'%2=(0 1 0) and so, q2 = yo. 
3. y3 = £3 — 191.43 = (0 0 1)° andso, g3 = y3. 


We have atrived at 


1 

0 1 
. Once the idea is grasped the actual calculations are best left to a machine. 
Matlab accomplishes this via the orth command. Its implementation is a bit 


more sophisticated than a blind run through our steps (1) through (m). As a 
result, there is no guarantee that it will return the same basis. For example 


>>X=[1 11;011;0 0 1]; 
>>Q=orth(X) 
Q= 

0.7370 -0.5910 0.3280 
0.5910 0.3280 -0.7370 


0.3280 0.7370 0.5910 


This ambiguity does not bother us, for one orthogonal basis is as good as 
another. Let us put this into practice, via (10.8). 


The Diagonalization of a Symmetric Matrix 


By choosing an orthogonal basis {q;,4.|1 << k < n;} for each R(P;) and 
collecting the basis vectors in 


Qj = (Gr 9,2 +++ In;) 
we find that 
nj 
P,=0)0;7 =) anaje’ 
k=1 
As a result, the spectral representation takes the form 
Equation: 
h T 
B= 37j-1459jQ; 
h nN; 
= D4 Aj pe Vik 


This is the spectral representation in perhaps its most detailed dress. There 
exists, however, still another form! It is a form that you are likely to see in 
future engineering courses and is achieved by assembling the Q; into a 
single n-by-n orthonormal matrix 


QO=(Oi ax :Op) 


Having orthonormal columns it follows that Q7Q = I. Q being square, it 
follows in addition that Q7 = Q-!. Now 


Baik = 255k 


may be encoded in matrix terms via 
Equation: 


BQ=QA 


where A is the n-by- n diagonal matrix whose first n; diagonal terms are 
Ai, whose next nz diagonal terms are Ag, and so on. That is, each 4; is 
repeated according to its multiplicity. Multiplying each side of [link], from 
the right, by Q7 we arrive at 

Equation: 


B=QAQ? 


Because one may just as easily write 
Equation: 


QTBQ=A 


one says that Q diagonalizes B. 


Let us return the our example 
111 
B= iid 
tf, 1 a 


of the last chapter. Recall that the eigenspace associated with A, = 0 had 


=—1 
ai= 1 


and 


—l 
€12 = 
1 


for a basis. Via Gram-Schmidt we may replace this with 


: 1 
i 
ve 4 
and 
—1 
: 1 
Tio Sa 
V6 5 
Normalizing the vector associated with Az = 3 we arrive at 
i 1 
Qi=— | 
Vo 5 
and hence 
-J/3 -1 v2 
1 
Q= (4 q) @) == V3 -1 V2 
6 
0 @ we 
and 


The Matrix Exponential 


The matrix exponential is a powerful means for representing the solution to 7n linear, 
constant coefficient, differential equations. The initial value problem for such a system 
may be written 


a'(t) = Az(t) 
x(0) = x 


where A is the n-by-n matrix of coefficients. By analogy to the 1-by-1 case we might 
expect 
Equation: 


a(t) = e“’u 


to hold. Our expectations are granted if we properly define e4*. Do you see why simply 
exponentiating each element of At will no suffice? 


There are at least 4 distinct (but of course equivalent) approaches to properly defining 
e4t. The first two are natural analogs of the single variable case while the latter two 
make use of heavier matrix algebra machinery. 


1. The Matrix Exponential as a Limit of Powers 

2. The Matrix Exponential as a Sum of Powers 

3. The Matrix Exponential via the Laplace Transform 

4. The Matrix Exponential via Eigenvalues and Eigenvectors 


Please visit each of these modules to see the definition and a number of examples. 


For a concrete application of these methods to a real dynamical system, please visit the 
Mass-Spring-Damper module. 


Regardless of the approach, the matrix exponential may be shown to obey the 3 lovely 
properties 


1; 4 (e“*) Ac Se" A 
2. eAltitty) — ettieAh 


3. e“¢ is nonsingular and (e4t)* = e (ad) 


Let us confirm each of these on the suite of examples used in the submodules. 


Example: 


If 
I 0 
A= 
(0 2) 
then 
0) 
“6 3) 
0 eZ! 
s 1 0 ee 8 
Chalo eee = : 
QO 2e* QE2 Que 
elitte 0 _ ete” 0 7 Ca) e2? 0 
2. ft) e2tit 2te a ) e2tie2te =, ) e2ti ) e2t2 
—t 
Ay Ve eed 0 ay (CAP 
3. (e ) =((, ye ) 
Example: 
If 
GQ 1 
A= 
Ga) 
then 


—sin(t) cos(t 


. 4 (eAt) _ — sin(t) os ) ye oe cole) ) 
é —cos(t) —sin(t) —cos(t) —sin(t) 
2. You will recognize this statement as a basic trig identity 
( cos(t, + ta) ee = ( cos(t;) ao) ( cos(tg) on 


—sin(t, + tg) cos(t, + t) —sin(t,) cos(t,;)/ \—sin(t2) cos(t) 


wee cos(t) ae ) 


3. Using the cofactor expansion we find 


os ee a) = (a t) 


sin(—t) 


sin(t) cos(t) 


Example: 
If 
then 


1u+e\ (fl i\ fl t& 

0 1 TN melee Ore a 
—l1 

& ‘) =((; ieee 

oy al Cea 


sin(—t) 
cos(—t) 


) 


= e (At) 


The Matrix Exponential as a Limit of Powers 
You may recall from Calculus that for any numbers a and t one may 


achieve e™ via 
Equation: 


The natural matrix definition is therefore 
Equation: 


where J is the n-by-n identity matrix. 


Example: 
The easiest case is the diagonal case, e.g., 


a=(03) 


for then 


Note that this is NOT the exponential of each element of A. 


Example: 
As a concrete example let us suppose 


From 
Ag 
I+ At= 
—t l 
(x ) 1 £f i¢ i 
7a aes, =4 = 2 
2 ze oroee = 1-45 
3 t 3 
(1+) = 1-3 t-y 
t? t? 
3 -t+£ 1-4 
At\'’ 4+ +1 _£ 
(1+) = 8 256 16 
as 
: itt bog tl 
—— ie 26 t 
cae) = + a95 +1 — 95 1 3125 
243 t? —- 
y eee toe 


We discern a pattern: the diagonal elements are equal even polynomials 
while the off diagonal elements are equal but opposite odd polynomials. 
The degree of the polynomial will grow with k and in the limit we 
recognize’ 


ie ( cos(t) — 


—sin(t) cos(t) 


Example: 
If 


then 


The Matrix Exponential as a Sum of Powers 


You may recall from Calculus that for any numbers a and t one may achieve 
e™ via 


Equation: 
k 
at ___ (at) 
= » kl 
k=0 
The natural matrix definition is therefore 
Equation: 
At ow (At)" 
=) kl 


where A® = J is the n-by-n identity matrix. 


Example: 
The easiest case is the diagonal case, e.g., 


a=(93) 


Zor = 


for then 


and so (recalling [link] above) 


Note that this is NOT the exponential of each element of A. 


Example: 


As a second example let us suppose 


(49) 


We recognize that its powers cycle, i.e., 


(peak 
a =A 
(Ao) 
and so 
2. 4 3 5 
oft 1-54+G-... t-Fta-... _ ( cos(t) — sin(t) 
Pg Se ee —sin(t) cos(t) 
Example: 
If 


The Matrix Exponential via The Laplace Transform 


You may recall from the Laplace Transform module that may achieve e” 
via 
Equation: 


The natural matrix definition is therefore 
Equation: 


et Y! (sI— A)" 


where J is the n-by-n identity matrix. 


Example: 
The easiest case is the diagonal case, e.g., 


“63 


for then 
1 
0 
—1 s—l 
(sf — A) = ; i 
s—2 
and so (recalling [link] above) 
Sint 
0 go = Qe 


Example: 
As a second example let us suppose 


and compute, in matlab, 


>> inv(s*eye(2)-A) 


ans = [ sS/(S42+1), 1/(sS42+1) ] 
[-1/(s42+1), s/(S2+1) ] 


>> ilaplace(ans) 


ans = [ cos(t), sin(t)] 
[-sin(t), cos(t)] 


Example: 
If 


a=(5 0) 


then 


>> 


>> 


inv(s*eye(2)-A) 


ans = [ 1/s, 1/s42] 
[ 0, 1/s] 


ilaplace(ans) 


ans = [ 1, t] 
Lo 4] 


The Matrix Exponential via Eigenvalues and Eigenvectors 


In this module we exploit the fact that the matrix exponential of a diagonal 
matrix is the diagonal matrix of element exponentials. In order to exploit it 
we need to recall that all matrices are almost diagonalizable. Let us begin 
with the clean case: if A is n-by-n and has n distinct eigenvalues, \;, and 
therefore n linear eigenvectors, s;, then we note that 


VIE {1, ee .,n} : (As, = Aj 8;) 


may be written 


Equation: 
A= SAS" 
where S = ($1 82 ... $n) is the full matrix of eigenvectors and 
A = diag (Ai, A2,..-, An) is the diagonal matrix of eigenvalues. One cool 


reason for writing A as in [link] is that 
AP aShs “Sis "= so 
and, more generally 
A= SA's 


If we now plug this into the definition in The Matrix Exponential as a Sum 
of Powers, we find 


eft -~ Se4tg-1 
h At - P 
where e*~ is simply 
diag(e™ 6" ane") 


Let us exercise this on our standard suite of examples. 


Example: 


If 
1 0 
A= 
(0 2) 

then 

SAA 
and so e4? = e“*. This was too easy! 
Example: 


As a second example let us suppose 


and compute, in matlab, 


>> [S, Lam] = e1ig(A) 


S = 0.7071 0.7071 
0 + 0.70711 0 - 0.70711 
Lam = 0 + 1.00001 0 
0 0 - 1.00001 


>> Si = inv(S) 


0.7071 0 - 0.70711 
0.7071 Oot On 7OFi4 


ok 


>> simple(S*diag(exp(diag(Lam)*t))*S1) 


ans = [ cos(t), Sin(t)] 
[-sin(t), cos(t)] 


Example: 
If 


then matlab delivers 


>> [S, Lam] = e1g(A) 


1.0000 -1.0000 
0 0.0000 


Lam = 0 0 


So zero is a double eigenvalue with but one eigenvector. Hence S is not 
invertible and we can not invoke ((link]). The generalization of ([link]) is 
often called the Jordan Canonical Form or the Spectral Representation. The 
latter reads 


h 
A Gee, 


j=l 


where the A, are the distinct eigenvalues of A while, in terms of the 
resolvent R(z) = (zI — A)™’, 


P; = — | R(z)d 
t 2mut Cg, ‘ 


is the associated eigen-projection and 


Be Record 


J Oni 


is the associated eigen-nilpotent. In each case, C’; is a small circle 
enclosing only Aj. 
Conversely we express the resolvent 


h mit 


HOS a 02 mae 


where 
m, —dim (2(P;)) 


with this preparation we recall Cauchy's integral formula for a smooth 
function f 


1 f f® 4 


ri zZ-a 


f(a) = 


where C(a) is a curve enclosing the point a. The natural matrix analog is 


A= = | HOR 


where C(r) encloses ALL of the eigenvalues of A. For f(z) = e” we find 
Equation: 


with regard to our example we find, h = 1, A; = 0, P; = I, m; = 2, 
D,; = Aso 


eee 


Let us consider a slightly bigger example, if 


aE 
As AL 
0 0 


wey oS S&S 


then 


>> R = inv(s*eye(3)-A) 


Pe a(p aber es ab sera. 0] 
L 0, WA SSI), 0] 
[ 0, 0, 1/(s-2)] 


and so A; = 1 and Ag = 2 while 


and som, = 2 


ol 0 0 
0 0 
and 
00 0 
Po Oy 00 
Oat 


and my = 1 and Dz = 0. Hence 


eft = e! (P; Se tD,) Sr e* P> 


The Mass-Spring-Damper System 
y, k 2c 
aa 


Mass, 
spring, 
damper 
system 


If one provides an initial displacement, xo, and velocity, vo, to the mass depicted in [link] then one 
finds that its displacement, x(t) at time t satisfies 


Equation: 
d? x(t) d x(t) 
ey P 1 ae T kax(t) = 0 
x0) = ao 
x'(0) = v9 


where prime denotes differentiation with respect to time. It is customary to write this single second 
order equation as a pair of first order equations. More precisely, we set 


ui(t) = x(t) 
w(t) = 2'() 


and note that [link] becomes 
Equation: 


Denoting u(t) = (w(t) u2(t))” we write [link] as 
Equation: 


We recall from The Matrix Exponential module that 


u(t) = e4tu(0) 


We shall proceed to compute the matrix exponential along the lines of The Matrix Exponential via 
Eigenvalues and Eignevectors module. To begin we record the resolvent 


ee —1 < i m ) 


mz2+2cze+k Mz 


The eigenvalues are the roots of mz” + 2cz + k, namely 


Vd,d= Vc?—mk: (» = —e) 


We naturally consider two cases, the first being 


1. d £ 0. In this case the partial fraction expansion of R(z) yields 


—-1 1 f{d-c -m —-1 1 /c+d m™ —1 —1 
(2) z= 3 ( k tos a (*, i.) as a ae 


At 


and so e4t = e*P, + e* Py. If we now suppose a negligible initial velocity, i.e., vo, it follows 
that 


Equation: 


x(t) = st (et (d—c) +e** (c+ d)) 


If dis real, i.e., if c? > mk, then both A; and Ag are negative real numbers and x(t) decays to 0 
without oscillation. If, on the contrary, d is imaginary, i.e., C< mk, then 
Equation: 


a(t) =e (cos(lat + i sn(\dt)) 


and so x decays to 0 in an oscillatory fashion. When [link] holds the system is said to be 
overdamped while when [link] governs then we speak of the system as underdamped. It 
remains to discuss the case of critical damping. 


.d = 0. In this case Ay = Ag = = and so we need only compute P; and Dy. As there is but 


N 


one P; and the P; are known to sum to the identity it follows that P,; = I. Similarly, this equation 
dictates that 


k 
4/ & 1 

D, = AP, —»1P, = A-—AT = - 
ee 


On substitution of this into this equation we find 


Equation: 


a (Vi) +t t 7 
- (+/£) 1—t4/£ 


Under the assumption, as above, that vg = 0, we deduce from [link] that 


x(t) = _ (v3) ( ia /*) x0 


The Singular Value Decomposition 


The singular value decomposition is another name for the spectral 
representation of a rectangular matrix. Of course if A is m-by-n and 

m # n then it does not make sense to speak of the eigenvalues of A. We 
may, however, rely on the previous section to give us relevant spectral 
representations of the two symmetric matrices 


© ATA 
e AAT 


That these two matrices together may indeed tell us ‘everything’ about A 
can be gleaned from 


N(A™A) = V(A) 

N (AAT) = W(AT) 
R(AtA) =R(A*) 
R(AA) = #A) 


You have proven the first of these in a previous exercise. The proof of the 
second is identical. The row and column space results follow from the first 
two via orthogonality. 


On the spectral side, we shall now see that the eigenvalues of AA? and 
A? A are nonnegative and that their nonzero eigenvalues coincide. Let us 
first confirm this on the adjacency matrix associated with the unstable 
swing (see figure in another module) 

Equation: 


The respective products are 


APA= 
= 


0 


oO fF & 
an) 
= © © © 


Analysis of the first is particularly simple. Its null space is clearly just the 
zero vector while A; = 2 and Az = 1 are its eigenvalues. Their geometric 
multiplicities aren; = 1 and ng = 2. In A’ Awe recognize the S matrix 
from the exercise in another moduleand recall that its eigenvalues are 

A, = 2, Ao = 1, and A3 = 0 with multiplicities ny = 1, ng = 2, and 

nz = 1. Hence, at least for this A, the eigenvalues of AA? and A? A are 
nonnegative and their nonzero eigenvalues coincide. In addition, the 
geometric multiplicities of the nonzero eigenvalues sum to 3, the rank of A. 
Proposition 


The eigenvalues of A.A? and A? A are nonnegative. Their nonzero 
eigenvalues, including geometric multiplicities, coincide. The geometric 
multiplicities of the nonzero eigenvalues sum to the rank of A. 


If A? Ax = Xx then x? A? Ar = dx" 7, ice., ( 
so A > 0. A similar argument works for AA’. 


| Aw ||)” = X(|| « ||)” and 


Now suppose that \; > 0 and that {z;,};’, constitutes an orthogonal 
basis for the eigenspace &(P;). Starting from 
Equation: 


T 
A Ax; k = AjL jk 


we find, on multiplying through (from the left) by A that 


AA" Ax; , = AJAX; % 


i.e., A; is an eigenvalue of AA” with eigenvector Azj;,, so long as 
Ax; # 0. It follows from the first paragraph of this proof that 

|| Ax; x ||= 4/2j, which, by hypothesis, is nonzero. Hence, 
Equation: 


is a collection of unit eigenvectors of AA” associated with \;. Let us now 
show that these vectors are orthonormal for fixed 7. 


T 1 ry T 
YViYi,k = 5 tiiA AX jb = L 5 4X5, = 0 
7 


We have now demonstrated that if \; > 0 is an eigenvalue of A’ A of 
geometric multiplicity n;. Reversing the argument, i.e., generating 
eigenvectors of A? A from those of A.A? we find that the geometric 
multiplicities must indeed coincide. 


Regarding the rank statement, we discern from [link] that if A; > 0 then 
Line A (A™A). The union of these vectors indeed constitutes a basis for 
& (AA), for anything orthogonal to each of these x ; ;, necessarily lies in 
the eigenspace corresponding to a zero eigenvalue, i.e., in VY (A™A). As 
RA’ A) = &(A") it follows that dim#(A* A) =r and hence the nj, 
for A; > 0, sum tor. 


Let us now gather together some of the separate pieces of the proof. For 
starters, we order the eigenvalues of A’ A from high to low, 


Ay > Ag >...>Ah 


and write 
Equation: 


ATA=XA,Xt 


where 
WX5 = {@jtye- er Vinyf ? (X= {X1,-- + Xa}) 


and A, is the n-by-n diagonal matrix with A, in the first n; slots, Az in the 
next 2 slots, etc. Similarly 
Equation: 


AA? =YApY* 


where 
VY; — {yjis - ar iigest : eg i= {¥1, 29 Yn}) 


and A,, is the m-by-m diagonal matrix with , in the first n; slots, Ag in 
the next 72 slots, etc. The y; , were defined in [link] under the assumption 
that A; > 0. If A; = 0 let Y; denote an orthonormal basis for. (AA‘). 


Finally, call 
2 oe 4/ A j 


and let +’ denote the m-by-n matrix diagonal matrix with a in the first n1 
slots, 72 in the next 7m» slots, etc. Notice that 
Equation: 

STS =A, 


Equation: 


Spee 


Now recognize that [link] may be written 


ALi k = O7Y5,k 
and that this is simply the column by column rendition of 
AX = YS; 


As X XT = I we may multiply through (from the right) by X7 and arrive 
at the singular value decomposition of A, 
Equation: 


A=YXX?t 


Let us confirm this on the A matrix in [link]. We have 


2000 
0100 
Age 

0010 

00 0 
4 © © 4 

: _ 
yo. V2 0 0 
V2 0 01 
0 v2 0 


Hence, 
Equation: 


andso A = YX'X? says that A should coincide with 


_ ov 

01 0\ /v2 000 
0 1 0 0 

100 0100 
000 1 
a ae 0 


This indeed agrees with A. It also agrees (up to sign changes on the 
columns of X) with what one receives upon typing [Y, SIG, X] = 
scd(A) in Matlab. 


You now ask what we get for our troubles. I express the first dividend as a 


proposition that looks to me like a quantitative version of the fundamental 
theorem of linear algebra. 
Proposition 


If YX" is the singular value decomposition of A then 


1. The rank of A, call it r, is the number of nonzero elements in ». 

2. The first r columns of X constitute an orthonormal basis for Z (A7). 
The n — rlast columns of X constitute an orthonormal basis for 
MN (A). 

3. The first r columns of Y constitute an orthonormal basis for #(A). 
The m — r last columns of Y constitute an orthonormal basis for 


WN (AP), 


Let us now 'solve' Aw = b with the help of the pseudo-inverse of A. You 
know the 'right' thing to do, namely reciprocate all of the nonzero singular 
values. Because m is not necessarily n we must also be careful with 
dimensions. To be precise, let /’* denote the n-by-m matrix whose first 4 


diagonal elements are 1 whose next n9 diagonal elements are = and so 


oy 
on. In the case that o;, = 0, set the final np, diagonal elements of X'* to 
zero. Now, one defines the pseudo-inverse of A to be 


At=xxytyt 


In the case of that A is that appearing in [link] we find 


J/2 00 
yt+_| 9 1 0 
0 01 
0 0 0 
and so 
=1 1 
= 0 ~ 0 ak. 
ve; we ve" °\ 19 1 9 
0 1 0 0 
At = Ye) Ea OO 
0 0 0 1 0 01 001 
1 90 + 0 0 00 
V2 V2 
therefore, 
-1 
0 = 0 
1 0 O 
At= ; 
0 5 O 
0 0 1 


in agreement with what appears from pinv(A). Let us now investigate the 
sense in which A‘ is the inverse of A. Suppose that b € R™ and that we 
wish to solve Ax = b. We suspect that A*b should be a good candidate. 
Observe by [link] that 


(ATA) AbD = XA,X*°XITY"D 


because X7X = J 
(AT A)Atb= XA, D*Y*D 
by [link] and [link] 
(ATA) Atb= XZ7ITLtY*D 
because S27 Ut = YF 
(ATA) Atb= XE" Y"D 
by [link] 
(ATA) Atb = AT 


that is A*b satisfies the least-squares problem A? Ax = A‘ b. 


